Electronic properties of bilayer and multilayer graphene 
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We study the effects of site dilution disorder on the electronic properties in graphene multilayers, 
in particular the bilayer and the infinite stack. The simplicity of the model allows for an easy 
implementation of the coherent potential approximation and some analytical results. Within the 
model we compute the self-energies, the density of states and the spectral functions. Moreover, we 
obtain the frequency and temperature dependence of the conductivity as well as the DC conductivity. 
The c-axis response is unconventional in the sense that impurities increase the response for low 
enough doping. We also study the problem of impurities in the biased graphene bilayer. 
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I. INTRODUCTION 

The isolation of single layer graphene by Novoselov et 
al.^ has generated enormous interest in the physics com- 
munity On the one hand, the electronic excitations of 
graphene can be described by the two-dimensional (2D) 
Dirac equation, creating connections with certain theo- 
ries in particle physicSfS, Moreover, the "relativistic" na- 
ture of the quasiparticles, albeit with a speed of prop- 
agation, vp, 300 times smaller than the speed of light, 
leads to unusual spectroscopic, transport, and thermo- 
dynamic properties that are at odds with the standard 
Landau- Fermi liquid theory of metals 4 On the other 
hand, graphene opens the doors for an all-carbon based 
micro-electronics,^ 

Due to the strong nature of the a bonds in graphene, 
and strong mechanical stability of the graphene lattice, 
miniaturization can be obtained at sizes of order of a few 
nanometers, beyond what can obtained with the current 
silicon technology (the smallest size being of the order 
of the benzene molecule). Furthermore, the same sta- 
bility allows for creation of entire devices (transistors, 
wires, and contacts) carved out of the same graphene 
sheet, reducing tremendously the energy loss, and hence 
heating, created by contacts between different materials.^ 
Early proposals for the control of the electronic prop- 
erties in graphene, such as the opening of gaps, were 
based on controlling its geometry, either by reducing it 
to nanoribbons,^ or producing graphene quantum dotst^ 
Nevertheless, current lithographic techniques that can 
produce such nanostructures do not have enough accu- 
racy to cut graphene to Angstrom precision. As a result, 
graphene nanostructures unavoidably have rough edges 
which have strong effects in the transport properties of 
nanoribbons.* In addition, the small size of these struc- 
tures can lead to strong enhancement of the Coulomb 
interaction between electrons which, allied to the dis- 
order at the edge of the nanostructures, can lead to 
Coulomb blockade effects easily observable in transport 



and spectroscopy^ 

Hence, the control of electronic gaps by finite geome- 
try is still very unreliable at this point in time and one 
should look for control in bulk systems which are in- 
sensitive to edge disorder. Fortunately, graphene is an 
extremely flexible material from the electronic point of 
view and electronic gaps can be controlled. This can 
be accomplished in a graphene bilayer with an electric 
field applied perpendicular to the plane. It was shown 
theoreticallyiSiii and demonstrated experimentall y^^i^^ 
that a graphene bilayer is the only material with semi- 
conducting properties that can be controlled by electric 
field effect. The size of the gap between conduction and 
valence bands is proportional to the voltage drop between 
the two graphene planes and can be as large as 0.1 — 0.3 
eV, allowing for novel terahertz devices^^ and carbon- 
based quantum dots^^ and transistors.— 

Nevertheless, just as single layer graphene ,-12. bilayer 
graphene is also sensitive to the unavoidable disorder 
generated by the environment of the Si02 substrate: 
adatoms, ionized impurities, etc. Disorder generates a 
scattering rate r and hence a characteristic energy scale 
h/r which is the order of the Fermi energy Ep = Tivpkp 
(kp cx is the Fermi momentum and n is the pla- 
nar density of electrons) when the chemical potential is 
close to the Dirac point [n 0). Thus, one expects 
disorder to have a strong effect in the physical proper- 
ties of graphene. Indeed, theoretical studies of the ef- 
fect of disorder in unbiasedii and biased'^^ graphene bi- 
layer (and multilayer) show that disorder leads to strong 
modifications of its transport and spectroscopic prop- 
erties. The understanding of the effects of disorder in 
this new class of materials is fundamental for any future 
technological applications. In this context it is worth 
to mention the transport theories based on the Boltz- 
mann equatioufi^i^S a study of weak localization in bi- 
layer graphene,— and also corresponding further experi- 
mental characterization] ^^1^'^ DC transport in few-layer 
graphene systems have been studied in Ref. [U both 
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without and in the presence of a magnetic field. 

In this paper, we study the efi^ects of site dilution 
(or unitary scattering) on the electronic properties of 
graphene multilayers within the well-known coherent po- 
tential approximation (CPA). While the CPA does not 
take into account electron localization , ^^'^^ it does pro- 
vide quantitative and qualitative information on the ef- 
fect of disorder in the electronic excitations. Further- 
more, this approximation allows for analytical results 
of electronic self-energies, allowing us to compute phys- 
ical quantities such as spectral functions (measurable 
by angle resolved photoemission, ARPES— i^^i^^i^i^) 
and density of states (measurable by scanning tunnel- 
ing microscopy, STM^i^i^ii^) , besides standard trans- 
port properties such as the DC and AC conductivities 
Furthermore, in the case of the semi-infinite stack of 
graphene planes we can compute the c-axis response of 
the system which is rather unusual since it increases with 
disorder at low electronic densities, in agreement with 
early transport measurements in graphite.— 

The paper is organized as follows. In Sec. |TT] we dis- 
cuss the band model of the graphene bilayer within the 
tight-binding approximation. We also connect our no- 
tation with the one established for graphite, namely 
the Slonsczewki-Weiss-McClure (SWM) parameteriza- 
tion. In Sec. mil we introduce several simplified band 
models and compare the electronic bands in different ap- 
proximations. The Green's functions that we use later 
on in the paper are given in Sec. IIVI 

We employ a simplified model for the disordered 
graphene bilayer in Sec. |V] and work out the conse- 
quences on the single particle properties encoded in the 
self- energies, the density of states (DOS) and the spec- 
tral function. Sec. IVII contains results for the graphene 
multilayer. In Sec. IVIll we introduce the linear response 
formulas that we use to calculate the electronic and op- 
tical response. The results for the conductivities in the 
bilayer are presented in Sec. IVIIH while those for the 
multilayer can be found in Sec. IIXI 

The rest of the paper concerns the problem of impu- 
rities in the biased graphene bilayer. The model of the 
system and some of its basic properties are discussed in 
Sec. |Xl In Sec. IXII we solve the problem of a Dirac delta 
impurity exactly within the effective mass approxima- 
tion. A simple estimate of when the interactions among 
impurities becomes important is presented in Sec. IXIII 
We treat more general impurity potentials with varia- 
tional methods in Sec. IXIIII and the special case of a 
potential well with finite range is studied in Sec lXIVl In 
Sec. IXVI we study the problem of a finite density of im- 
purities in the coherent potential approximation (CPA). 
The effects of trigonal distortions on our results for the 
biased graphene bilayer arc discussed briefly in Sec. lXVll 
Finally, the conclusions of the paper are to be found 
in Sec. IXVIII We have also included four appendices 
with technical details of the calculations of the mini- 
mal conductivity in bilayer graphene (App. the DOS 
in multilayer graphene (App. [B|, the conductivity ker- 



nels (App. [C]), and the Green's function in the biased 
graphene bilayer fApp. [D]). 



II. ELECTRONIC BANDS OF THE GRAPHENE 
BILAYER 

Many of the special properties of the graphene bilayer 
have their origin in its lattice structure that leads to the 
peculiar band structure that we discuss in detail in this 
section. A simple way of arriving at the band structure of 
the graphene bilayer is to use a tight-binding approxima- 
tion. The positions of the different atoms in the graphene 
bilayer are shown in Fig. [1] together with our labeling 
convention. 

The advantage of this notation is that one can discuss 
collectively about the A (B) atoms that are equivalent in 
their physical properties such as the weight of the wave 
functions and the distribution of the density of states etc. 
This notation was used in early work on graphitci^i^ 
Many authors use instead a notation similar to Al A, 
Bl B, A2 ^ B, and B2 A. In this notation the 
relative orientation within the planes of the A (A) and 
B (B) atoms are the same; but for the other physical 
properties the equivalent atoms are instead A (B) and B 
(A) . Because the other physical properties are often more 
relevant for the physics than the relative orientation of 
the atoms within the planes we choose to use the, in our 
view, most "natural" labeling convention. 



A. Monolayer graphene 

Let us briefly review the tight-binding model of mono- 
layer graphene*^ The band structure can be described 
in terms of a triangular lattice with two atoms per unit 
cell. The real-space lattice vectors can be taken to be 
ai = f(3, V3) andaa = f (3,-^3). Here a (« 1.4 A) de- 
notes the nearest neighbor carbon distance. Three vec- 
tors that connect atoms that are nearest neighbors are 
^1 = f(l,\/3), 52 = f(l,-V3), and ^3 = a(-l,0); we 
take these to connect the Al atoms to the Bl atoms. In 
terms of the operators that creates (annihilates) an elec- 
tron on the lattice site at position and lattice site aj 
[ a = (A, B) denotes the atom sublattice and j (j — 1) 
denotes the plane ] : c^^- ^, [c^- ^. ) , the tight-binding 
Hamiltonian reads: 

Ht.b.=tY^ (cai,r.Cbi3,+5, (1) 

Ri j = l,2,3 

Here t (w 3 eV) is the energy associated with the hopping 
of electrons between neighboring tt orbitals. We define 
the Fourier-transformed operators, 

Caj,R, ^ ^'^aj,k^ (2) 
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where A'" is the number of unit cells in the system. 
Throughout this paper we use units such that h = = 1 
unless specified otherwise. 

Because of the sublattice structure it is often con- 
venient to describe the system in terms of a spinor: 
in which case the Hamiltonian can 



be written as: 



t.b. 



C(k)A 

klC*(k) i 



(3) 



where 



C(k) = i^e' 



k<5. 



= te''^--/^ [2 cos(^^^^) + e-^3/c.a/2] 



(4) 



The reciprocal lattice vectors can be taken to be bi = 
§f(l,\/3) and ba = if(l>-\/3) as is readily verified. 
The center of the Brillouin zone (BZ) is denoted by 
r, but for the low-energy properties one can expand 
close to the K point of the BZ, which has coordinates 
K = 5^(0,-1). One then finds C(K+p) = a = uppe*", 
where vp — 3ta/2 and a = — eLTCteiii{px / Py) ■ Note that 
a = along the K — F line of the BZ and that it in- 
creases anti-clockwise. With these approximation one 
finds that the spectrum of Eq. ^ is that of massless 2D 
Dirac fermions: E± — ±v-pp. 




FIG. 1: [color online] Lattice structure of the graphene bi- 
layer. The A (B) sublattices are indicated by the darker 
(lighter) spheres and the planes are labeled by 1 and 2. 



B. Bilayer graphene 

Since the system is 2D only the relative position of 
the atoms projected on to the a;-y-plane enters into the 
model. The projected position of the different atoms are 
shown in Fig. O together with the BZ. Since the A atoms 
are sitting right on top of each other in the lattice, the 
hopping term between the Al and A2 atoms are local in 
real space and hence a constant that we denote by t± in 
momentum space. Referring back to Section FlI Al we note 



O Bl 




FIG. 2: [color online] The real space lattice structure of the 
graphene bilayer projected onto the x-y plane showing the 
relative positions of the different sublattices. The upper right 
corner shows the BZ of the graphene bilayer including the 
labeling of the high symmetry points. 



that the hopping Bl Al [Al Bl] gives rise to the 
factor C(k) [C*(k)], with C(k) defined in Eq. (g]). Since the 
geometrical role of the A and B atoms are interchanged 
between plane 1 and plane 2 we immediately find that in 
Fourier space the hopping A2 B2 [B2 A2] gives rise 
to the factor ((k) [C*(k)]. Furthermore, the direction in 
the hopping Bl B2 (projected on to the x-y plane) is 
opposite to that of hopping Bl Al. Thus we associate 
a factor W3C*(k) to the hopping Bl B2, where the 
factor W3 = 73/70 is needed because the hopping energy 
is 73 instead of 70 = t. Similarly, the direction of hopping 
Bl A2 (projected on to the x-y plane) is the same as 
Bl Al and therefore the term — U4C(k) goes with the 
hopping Bl — > A2. The minus sign in front of V4 follows 
from the conventional definition of 74 in graphite, as are 
discussed below. Continuing to fill in all the entries of the 
matrix the full tight-binding Hamiltonian in the graphene 
bilayer becomes: 



Ht.b.(k) = 



c 



c 


t± 


-V4C\ 


V/2 


-V4C 


f] 


-V4C 


-~V/2 + A 




V3C 


C 


-v/2/ 



(5) 



where the spinor is 5*1 = 



31. k' ''A2,k' '^B2,k) 



-k ^ ('^Al,k' "^E 

Here we have also introduced the conventional (from 
graphite) A that parametrizes the difference in energy 
between A and B atoms. In addition we included the 
parameter V which gives different values of the potential 
energy in the two planes, such a term is generally allowed 
by symmetry and is generated by an electric field that is 
perpendicular to the two layers. The system with V 
is called the biased graphene bilayer and has a gap in the 
spectrum, in contrast the spectrum is gapless if y = 0. 

It is also possible to include further hoppings into the 
tight-binding picture, this was done for graphite by John- 
son and Dresselhaus.^" The inclusion of such terms is nec- 
essary if one wants an accurate description of the bands 



4 



throughout the whole BZ. If we expand the expression 
in Eq. (O close to the K point in the BZ we obtain the 
matrix: 



where 



noip) 



V/2 + A 


a 






a* 


V/2 






t± 


—Via 


-V/2 + A 


a* 






a 


-V/2 



(6) 



where a was introduced after Eq. 

The typical behavior of the bands obtained from 
Eq. ^ is shown in Fig. [31 Two of the bands are moved 
away from the Dirac point by an energy that is ap- 
proximately given by the interplane hopping term t± for 

V ^ t±. In the figure we have taken V ^ 0; but for 

V = there is no gap for the two bands closest to zero 
energy (i.e. the Dirac point). 






p(eV) 



FIG. 3: [color online] Band dispersions near the K points in 
the bilayer along the direction a = 0, with V = 50meV and 
up = 1. Left: bands obtained from the full model in Eq. (|6]) 
with t± = 0.35 eV, = 0.1 and V4, — 0.05; Right: bands 
obtained from the simplified model in Eq. (|10[) . 



C. The SIonsczewki-Weiss-McClure (SWM) model 



First we make the observation that the graphene bi- 
layer in the A-B stacking is just the unit cell of graphite 
that we depict in Fig. [1] Therefore, if the two planes are 
equivalent much of the symmetry analysis of graphite is 
also valid for the graphene bilayer. Thus we could al- 
ternatively use the SWM for graphite with the proper 
identification of the parameters. The SWM model for 
graphite is usually written as 



"WswMC — 



El 





Hi3 


Hl3 \ 





E2 


EI23 


^H^3 


^13 


^23 


E3 


H33 


^^13 


~H23 


H33 


E3 J 



(7) 



El = A + 7ir + -75r2, 

E2 = A-7ir + i75r2, 

E3 = ^72^^ 

Hi3 = -^(-7o+74r)e*"C, 

H23 = -^(7o+74r)e"C, 

H33 = 73re'"C- 



(8a) 
(8b) 
(8c) 
(8d) 

(8e) 
(8f) 



Here ( — 3afc/2, and F = 2cos{k±d), with d « 3.7A be- 
ing the interplane distance. Typical values of the param- 
eters from the graphite literature are shown in Table U 



70 


71 


72 


73 


74 


75 


76 = A 


ep 


3.16 


0.39 


-0.02 


0.315 


0.044 


0.038 


0.008 


-.024 


3.12 


0.377 


-0.020 


0.29 


0.120 


0.0125 


0.004 


-.0206 



TABLE I: Values of the SWM parameters for the band struc- 
ture of graphite. Upper row from Ref. 36] and lower row from 
Ref. 0. 

It is straightforward to show that by identifying 71 = 
t± and taking 72 = 75 = 0, F = 1 and V = the 
matrices in Eq. Q and Eq. ^ are equivalent up to a 
unitary transformation. Hence they give rise to iden- 
tical eigenvalues and band structures. This completes 
the correspondence between the tight-binding model and 
the SWM model (see also Refs. [iOld^l for a discussion 
on the connection between the tight-binding parameters 
and those of SWM.) 

The accepted parameters from the graphite literature 
results in electrons near the K point [k± — 0] and holes 
near the H point [k± = 7r/(2c?)] in the BZ as sketched 
in Fig. m These electron and hole pockets are mainly 
generated by the coupling 72 that in the tight-binding 
model corresponds to a hopping between the B-atoms of 
next-nearest planes. Note that this process involves a 
hopping of a distance as large as ~ 7 A. 

Finally, it is interesting to note that at the H-point in 
the BZ, F = 0, and therefore the two planes "decouple" at 
this point. Furthermore, if one neglects A the spectrum 
is that of massless Dirac fermions just like in the case 
of graphene. Note that in graphite A and B atoms are 
different however, and that the term parametrized by A, 
that breaks sublattice symmetry in each plane, opens a 
gap in the spectrum leading to massive Dirac fermions 
at the H-point. Since the value of A in the literature is 
quite small, the almost linear massless behavior should 
be observed by experimental probes that are not sensitive 
to this small energy scale. 

The values of the parameters used in the graphite liter- 
ature are consistent with a large number of experiments. 
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FIG. 4: [color online] Left: graphite lattice; Right: three- 
dimensional BZ with the symmetry points K and H indicated. 
The accepted parameters for graphite results in electron pock- 
ets near the K points and hole pockets near the H points as 
sketched in the figure. 



The most accurate ones are various magneto-transport 
measurements discussed in Ref. [s^. More recently, 
angle-resolved photo-emission spectroscopy (ARPES) 
was used to directly visualize the dispersion of massless 
Dirac quasi-particles near the H point and massive quasi- 
particles near the K point in the BZ i^'^'^^i'^° 

The band structure of graphite has been calculated and 
recalculated many times over the years, a recent reference 
is Ref. [1^. It is also worth to mention that because of 
the (relatively) large contribution of the non-local van 
der Waals interaction to the interaction between the lay- 
ers in graphite, the usual local density approximation 
or semilocal density approximation schemes are off by 
an order of magnitude when the binding energy of the 
planes are calculated and compared with experiments. 
For a discussion of this topic and a possible remedy, see 
Ref. El. 



III. SIMPLIFIED ELECTRONIC BAND 
MODELS 

In this section we introduce three simplified models 
that we employ for most of the calculations in this paper. 
We also show how to obtain an effective two-band model 
that is sometimes useful. 



near the K point in the BZ. Here we write + iky = 

^g«0(k) ^ -vvhere k = yk^ + k^ and (/>(k) is the appropriate 

angle. Note that the absolute value of the angle can 
be changed by a gauge transformation into a phase of 
the wave functions on the B sublattices. This reflects 
the rotational symmetry of the model. If one includes 
the "trigonal distortion" term parametrized by 73 the 
rotational symmetry is broken and it is necessary to keep 
track of the absolute value of the angle. From now on in 
this paper, we most often use units such that vp = 1 for 
simplicity. 

This Hamiltonian has the advantage that it allows for 
relatively simple calculations. Some of the fine details of 
the physics might not be accurate but it works as a mini- 
mal model and capture most of the important physics. It 
is important to know the qualitative nature of the terms 
that are neglected in this approximation, this is discussed 
later in this section. It is also an interesting toy model as 
it allows for (approximately) "chiral" particles with mass 
(i.e., a parabolic spectrum) at low energiesJ^ 



B. Biased graphene bilayer 

For the biased bilayer, a minimal model employs 
Eq. ([9]) augmented with the bias potential V: 



( 



\ 



VII 




(k) 



F/2 







ke 



k) 



-V/2 J 



(10) 

This model was introduced in Refs. |llll45l |. It correctly 
captures the formation of an electronic gap of size ~ V 
at the K point and the overall features of the bands as 
can be seen in Fig. [3] Nevertheless, the fine details of the 
bands close to the band edge are not captured correctly 
in this simple model; this fact is illustrated in Figs. [S][S1 
In particular the simple model is cylindrically symmetric; 
whereas the "trigonal distortion" breaks this symmetry. 
Thus the inclusion of V3 results in a "trihorn" structure 
for small values of V and a weaker modulation of the 
band edge for larger values of V as illustrated in Fig. [51 



A. Unbiased bilayer 



C. Multilayer graphene 



For the unbiased bilayer, a minimal model includes 
only the nearest neighbor hopping energies within the 
planes and the intcrplanc hopping term between A 
atoms; this leads to a Hamiltonian matrix of the form: 



f 



HB(k) 









ke'"!''^^ 

ke"^^^'> 



(9) 



In the graphene multilayer, a minimal model for the 
bands is again given by Eq. ^ with the understanding 
that the momentum label also includes the perpendicular 
direction: k — > (ky ,kj_). The only change is that we must 
make the substitution t± ~-^ 2t± cos{k±d) everywhere. 
Note that this is exactly the F-factor appearing in the 
SWM model discussed in Sect ion Hi CI In the following we 
often use units such that the interplane distance d is set to 
1, then ~ since the unit cell holds two layers - the allowed 
values of k± lies in the interval [— 7r/2, 7r/2]. We note 
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FIG. 5: [color online] Band dispersions near the band edge 
(note the energy scale) in the biased graphene bilayer. Left: 
V = 50meV, Right: V = 200 meV. The solid line is the 
simplified model in Eq. (|10p that is cylindrically symmetric. 
The other lines are along difi^erent directions in the BZ for the 
full model in Eq. ©: a = (dotted), a = tt/9 (dash-dotted), 
and a — 2n/9 (dashed). The parameters are the same as in 
Fig. El 




FIG. 6: [color online] Contour plots of the band dispersions 
near the band edge in the biased graphene bilayer for V — 
50meV. Left: full model in Eq. ([6}, Right: simplified model 
in Eq. pO[) . The parameters are the same as in Fig. |3] 



that this band model was used already in the seminal 
paper by Wallace as a simple model for graphite. '^^ More 
recent works on the band structure of few-layer graphene 
systems include Refs. mHililil 



D. Approximate effective two- band models 

There are two main reasons for constructing approxi- 
mate two-band models. Firstly, on physical grounds the 
high-energy bands (far away from the Dirac point) should 
not be very important for the low-energy properties of 
the system. Secondly, it is sometimes easier to work with 
2x2 instead of 4 x 4 matrices. Nonetheless, it is not al- 
ways a simplification to use the 2-band model when one is 
studying inhomogeneous systems as it generally leads to 
two coupled second order differential equations whereas 
the 4-band model involves four coupled linear differential 
equations. The matching of the wave functions in the 2- 
band case then involves both the continuity of the wave 
function and its derivative whereas in the 4-band model 
only continuity of the wave function is necessary. We 



note that the two-band description of the problem is only 
valid as long the electronic density is low enough, that 
is, when the Fermi energy is much smaller than t±. At 
intermediate to high densities a 4-band model is required 
in order to obtain the correct physical properties'^. 

In this section, we derive the low-energy effective model 
by doing degenerate second order perturbation theory. 
The quality of the expansion is good as long as vpp <C 
t± w 0.35 eV. We first present the general expression for 
the second-order 2x2 effective Hamiltonian, thereafter 
various simplified forms are introduced. Analyses similar 
to the one here were previously described in Refs. Tolli^ . 

First we introduce the projector matrices Vq ~ 
Diag[0, 1, 0, 1] (Pi = Diag[l, 0, 1, 0]) that projects onto 
(out of) the low-energy subspace of the B atoms. Then 
we split the Hamiltonian in Eq. ([6]) according to: Hq = 
K-Q + K-i + IC2, with 



/Co 



/Ci = vf 



fA + V/2 







\ 







V/2 












t± 


A 


- v/2 









\ 








~V/2) 




/ 


pe"t' 







-Vipe 







—V4pe 












-v^pe"*' 







pe-' 


\-Vipe"'' 





pe^'f' 











/o 





^ 











U3Pe"^ 





















\0 W3pe~' 


4> 


) 







(11) 



(12) 



(13) 



Introducing the vectors {10} that to zeroth order only 
have components in the low energy subspace (i.e., 
Vi\l^°^) = 0) and following the standard procedure (see 
e.g., Ref. [sd j) for degenerate perturbation theory we ar- 
rive at: 



{i\ronoPo\i') ~ {i\Vo[iCo + x^iC2]Vo\i') 

1 



E-ICo 



VitCiVoll'), (14) 



where we explicitly assume that IC2 is of the same order 
as K-l/K-Q. This expression is correct to second order in A. 
Note that we are doing second order perturbation theory 
for all of the components of the 2 x 2-matrix in the low- 
energy subspace. Working to first order in (and then 
setting A = 1) one obtains for this matrix (taking vp = 1 
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for brevity): 



^low 



y/2 V3pe"l' 



— i<p 



-V+2t^Vi+A{l+vi) 
t^-A(A-y) 
t^(l+i'4)+2f4A 2»<j) 



V+2txf4+A(l+ii^) 
tl-A(A+y) 



(15) 



Taking A = leads to an even simpler expression, in 
particular the effective spectrum becomes: 



£^low,± 



± 



2p2 1 2 y2 



2v3P^[t^{l+vl) 



cos(30). 



(16) 



r 



That this approximation to the bands is excellent near 
the band edge for small values of the bias V is illustrated 
in Fig. [71 For larger values of the bias the agreement 
is less accurate because the assumption of smallness of 
certain terms in the perturbation expansion is no longer 
valid. 
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It is also easy to co mpute the energy eigenval ues that are 
given by tj2 ± v'ii/4 + fc2 and ± ^iy4 + fc2. 

Before we solve for the Green's functions it is convenient 
to allow for a local frequency dependent self-energy in 
the problem. In the general case the self-energies on all 
of the inequivalent sites in the problem are allowed to be 
different, and we explicitly introduce the matrix 



(18) 



p(eV) 



\\\ If 

'■V \i 




/Sai(cc') \ 


'.\\ '/■'/ 




SBi(t^) 

SA2(t^) 






^ ^B2{UJ)I 




to describe this. The Green's function matrix 


-0.2 -0.1 0.1 0.2 

p(eV) 


given by the equation 



FIG. 7: [color online] Band dispersions near the band edge in 
the biased graphene bilayer along two directions in the BZ. 
Left: V = 50meV, Right: V = 200 meV. The solid (dash- 
dotted) line is the effective model in Eq. ((16} along a = 
(a = 7r/6). The dotted (dashed) line is the full model in 
Eq. (|6]) along a = (a = 7r/6). The parameters are the same 
as in Fig. [S] For V = 50 meV the different curves are almost 
not discernible. 



G-i(cc;,k) =w-7io(k)-7is(c 



(19) 



In the case of the unbiased bilayer the A (B) sites in 
both of the layers are equivalent and we only need two 
self-energies: and S]b('-^) which are local but we 

allow for a frequency (w) dependence. In this case the 
matrix inversion is simple since it factorizes into two 2x2 
matrices. An explicit form is given by 



GREEN'S FUNCTION IN THE GRAPHENE 
BILAYER 



IV. 



As discussed in Section Hill we use the minimal model 
Hamiltonian in Eq. We note that the phases (j) = 
(^(k) can be gauged away by an application of a unitary 
transformation defined by the matrix 



A^i(k) 



/l \ 

e-'-^e^) 

10 

\0 e"f'^^^) 



(17) 



G(.,k) . A^,(k) Q^'S] ^^(^)' 

\^g^^^(w,fc) g"{uj,k) J 

(20) 

where k — |k|. Here D (ND) stands for diagonal (non- 
diagonal) in the layer index. The components of the g- 
matrices are given by 



D,ND 



D.ND 



D.ND 



2D^ 



2D. 



(21a) 



--t^-^A^^^^ + tL^^ (21b) 



2D- 
k ^ k 



2D4 



2D- 2D4 



(21c) 
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where 



(22) 



Note that we often suppress the momentum and fre- 
quency dependence in the foUowing when no confusion 
arises. We wih come back to the biased case in Sec- 
tion |XVl 



V. IMPURITY SCATTERING: T-MATRIX AND 
COHERENT POTENTIAL APPROXIMATION 

We are interested in the influence of disorder in the 
bilayer. To model the impurities we use the standard 
t-matrix approach and the Coherent Potential Approxi- 
mation (CPA). The effect of repeated scattering from a 
single impurity can be encoded in a self-energy which can 
be computed fromi^ 



VjNGVj + VjNGVjNGVj 



(23) 



Here Vj is a matrix that encodes both the strength and 
the lattice site of the impurity in question. For example, 
an impurity on an Al lattice site of strength U at the 
origin is encoded in Fourier space by the matrix 



V,= — 



/l o\ 



Vo oy 



(24) 



implying that the potential is located only on a single 
site. We have also introduced the quantity: 



27r 



G„,„j(a;,k), (25) 



which is the local propagator at the impurity site; and 
in the second step the k-sum is to be taken over the 
whole BZ. The last line is an approximate expression 
that is obtained by expanding the propagator close to 
the K points and taking the continuum limit with the 
introducing of the cutoff A. We estimate the cutoff by a 
Debye approximation that conserves the number of states 
in the BZ. Then A w 7eV and in units of the cutoff we 
have t± Ri 0.05 (taking t± « 0.35 eV). Due to the special 
form of the propagator and the impurity potential the 
self-energy we get from this is diagonal. The result for 
site dilution (or vacancies) is obtained by taking the limit 
L/ ^ oo, so that the electrons are not allowed to enter 
the site in question. We also introduce a finite density 
rii of impurities in the system. To leading order in the 



impurity concentration the equations for the self-energies 
then become: 



-Gb. 



(26a) 
(26b) 



The explicit form of the propagators in Eq. (PT|) makes 
it easy to compute the G's. The t-matrix result for the 
self-energies is obtained by using the bare propagators on 
the right hand side of Eq. ([26]) . In the CPA one assumes 
that the electrons are moving in an effective medium with 
recovered translational invariance which in this case is 
characterized by the local self-energies. To determine 
what the medium is, one must solve the equations self- 
consistently with the full propagators on the right hand 
side of Eq. ((26l) . Because of the simple form of the propa- 
gators this is a simple numerical task in the model we are 
using. To simplify the equations further we assume that 
A 3> w, Ea, Sb. This is a physical assumption since 
when the self-energies becomes of the order of the cut- 
off the effective theory breaks down. The self-consistent 
equations then reads: 



[^]-^ = -Ga 

Hi 



2A2 



J2 



a=± 



A2 



-{u; + atj^ - - Sb) 



, (27a) 



^ - Sa 
2A2 



A2 



-(cj - atj_ - Sa)(c^ - Sb) 
{u -tj_- ^a){(^ - Sb 



2A2 + t±- I]a)(cc - Sb) 



(27b) 



This includes inter-valley scattering in the intermediate 
states. It is easy to obtain the non-disordered density of 
states from these equations by taking Sa = Sb — and 
Lu ^ Lo + i6 (here (5 is a positive infinitesimal) resulting 



pUlo) ^ -i Im Ga = ^ [1 + &i\io\ - tj^)] , (28a) 



pUlu) ^ --ImGB 

TT 



(28b) 



Observe in particular that the density of states on the A 
sublattice goes to zero in the limit of zero frequency, this 
fact is responsible for much of the unconventional physics 
in the graphene bilayer. In contrast the density of states 
on the B sublattice is finite at w = 0. We discuss how 
this result is changed with disorder and the solution of 
Eq. (g?!) in the following. 
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A. Zero frequency limit 



self-energies is given by the contribution to second order: 



One interesting feature of the CPA equations in 
Eq. ((27l) is that it is easy to see that they do not al- 
low for a finite Sa in the limit of w ^ 0. Since by setting 
uj = the last term in Eq. (j27bp must vanish, and this 
is not possible for finite values of Sa- Then one also 
must have that Sb there. This implies that the 
density of states on sublattice A is still zero even within 
the CPA in the limit lu ^ 0. More explicitly, by defining 
Sa^b = — one can show (assuming Ea S> t± and 
Sb ^ that Sa and Sb are given asymptotically by: 



Sa 
Sb 



UiOJ J ' 



and ^ satisfies 



Clog(l/0- 



(29a) 
(29b) 

(30) 



Notice that \/^A ~ y^A is exactly the energy scale that 
is generated by disorder of the same kind in the single 
layer caseM^ We have checked that the expressions in 
Eq. (P^l seem to agree with the numerical calculations 
in the small frequency limit, and the frequency range in 
which they holds grows with increasing rii . 



B. Uncompensated impurity densities 

The divergence of the self-energy Sa on the A sub- 
lattice in the CPA in the above is due to the fact that 
there is a perfect compensation between the number of 
impurities on the two sublattices: rij^A = '^i,B- For the 
more general case where nj.A 7^ ni^B the divergence is not 
present so that Sa may become finite at w = 0. To make 
comparison with other work on the graphene bilayer it 
is fruitful to consider another extreme limit where only 
the B sites are affected by the disorder. Explicitly this 
means that we take ni^A — and B 7^ 0. The gener- 
alization of the CPA equations in Eq. ([^T]) for this case 
then immediately imply that TiAi^) = 0. In the limit of 
u! ^ 0, Sb is finite, purely imaginary, and given by: 



Sb(w = 0) 



2A2 



-iV. 



(31) 



C. Born scattering 

Another often studied limit is the one of weak im- 
purities, in particular Koshino and Ando have stud- 
ied electron transport in the graphene bilayer in this 
approximation.^^ This is the Born limit and it can be 
studied using perturbation theory in the strength of the 
impurities U . The leading non-trivial contribution to the 



Sa = UiUGAU, 
Sb = UiUGBU. 



(32a) 
(32b) 



If one substitutes the bare propagators on the right hand 
side one finds Sa = and Sb — —iTrnit±{U / A)"^ /2 at 
the Dirac point. Thus, to leading order Born scattering is 
formally equivalent to the previous case with vacancies on 
only sublattice B exactly at w = 0. The frequency range 
for which the uj = result is valid is different however. 



D. Self-energy comparisons and the density of 
states 

We compare the self-energies obtained from the t- 
matrix and the CPA. Within the t-matrix the self- 
energies are given by 



Sa = 
Sb = 



F0+Z7rp0(c.)' 
n. 



where the p^'s are given in Eq. 



and 



2A2 



log( 



A4 



FBO^-Re[G°B] = Fl + i^^lo. 



2A2 



(33a) 
(33b) 

, (34a) 
. (34b) 



The results for the self-energies in the two different ap- 
proximations are shown in Figs. [S] and O Note that at 
least on the scale of the figures the S^i diverges as w ^ 
in the CPA, as discussed above. The solution to the 
self-consistent equations also does not converge very well 
when they are pushed close to the limit of ijj — > 0. The to- 
tal DOS on the A-sublattice and B-sublattice is pictured 
in Fig.[Tni Note in particular that the case of rii = .0001 
closely resemble the non-interacting case except for the 
new low-energy feature. A possible interpretation of the 
enhancement of the DOS on the B sublattice near w = 
is in terms of the "half-localized" states (meaning they 
do not decay fast enough to be normalizable at infinity) 
that have been studied for monolayer graphene. Be- 
cause these states have weight on only one sublattice (the 
opposite one of the vacancy) the construction in Ref. [5^ 
is valid also in the graphene bilayer when there is a va- 
cancy on one of the A sublattices. For a discussion of 
the related problem of edge states in bilayer graphene 
see Ref. [HI. 



E. Spectral function 

The electron spectral function A(k,w), which is ob- 
servable in ARPES experiments, is defined by A(k, w) = 
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FIG. 8: [color online] Self-energies within the t-matrix ap- 
proximation in the bilayer as a function of the frequency to. 
Left: sublattice A; Right: sublattice B. 



FIG. 10: [color online] Local density of states p on the dif- 
ferent sublattices in the bilayer. Left: t-matrix; Right: CPA. 
Top: sublattice A; Bottom: sublattice B. 
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FIG. 9: [color online] Self-energies within the CPA in the 
bilayer as a function of the frequency lo. Left: sublattice A; 
Right: sublattice B. 



— Trace[Ini G(k, w)]/7r, so that in our case: 
A(k, c.) = - ^ {im [G^A (k, c.)] -f Im [GEb (k, c.)] } . (35) 

The spectral function in the fcxoj plane, calculated within 
the CPA, is pictured in Fig. [11] As can be seen in the fig- 
ures the low-energy branch becomes significantly blurred, 
especially for the higher impurity concentrations. Note 
also that the gap to the high-energy branch becomes 
slightly larger as the disorder value increases due to the 
fact that Sa is not negligible there. 

Examples of the momentum distribution curves 
(MDC's) and the energy distribution curves (EDC's) in 
the disordered graphene bilayer are shown in Figs. [12] 
and [T2] The evolution of the peaks from delta functions 
to broader peaks with increasing disorder is clear in the 
figure. 



= ii.tmui fi, = 0.0(1] 



FIG. 11: Intensity plot of the spectral function in the k x ui 
plane (normalized by the cutoff) in the bilayer for different 
impurity concentrations in the CPA approximation. From 
left to right: m = lO"*, 10"^ 5 x 10"^. 



VI. GREEN'S FUNCTION AND 
ONE-PARTICLE PROPERTIES IN MULTILAYER 
GRAPHENE 

We will use the extension of the bilayer model to the 
multilayer that we introduced in Section IIII CI As dis- 
cussed there we can immediately use the Hamiltonian in 
Eq. ^ with the understanding that the momentum label 
also includes the perpendicular direction: k — > (k||,fcj_) 
and by substituting t± 2t± cos{k±d) everywhere. In 
particular the Green's function including the self-energies 
are again given by the expressions in Eq. I|2ip and 
Eq. (|22)) with the substitution t± 2t± cos{k±d). 



A. Self-energies and the density of states 

To get the CPA equations we must evaluate the local 
propagator G that is given by the straightforward gener- 
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FIG. 12: MDC's in the graphene bilayer. The three panels 
are for diflerent impurity concentrations and are calculated in 
the CPA. From the top the energy cuts are at the energies: 
.0005, .0055, .0105, .0155, .0205, .0255, and .0305 in units of 
the cutoff A. The curves are uniformly displaced for clarity. 



From these equations one can easily obtain the non- 
interacting density of states by considering the clean re- 
tarded case and take Ea — T^b — and uj ^ uj + i6, from 
which we get: 



1^1 
7rA2 



-+tan \ 



) Qi2t^-\Lo\) 



+ Lie(|c.|-2i^), (38a) 



Pb ^ Pa 



7rA2 



A2 



e{2ti_ - \uj\ 



(38b) 



The equivalent expression were previously obtained in 
Ref. [45| using a different method. The self-energy within 
the t-matrix is again given by the expression in Eq. (j33p , 
with the p's given by the non-interacting density of states 
in Eq. ([55]) . The F's are obtained from the real parts of 
the non-interacting local propagators: 
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FIG. 13: EDO's in the graphene bilayer. The three panels 
are for different impurity concentrations and are calculated 
in the CPA. From the top the cuts are at fixed values of the 
in-plane momentum k: .0001, .01, .02, .03, .04, .05, and .06 
in units of the cutoff A. The curves are uniformly displaced 
for clarity. 



alization of Eq. (|25p to include an extra sum over k± : 

GH = ^Ei^EGM- (36) 



Here Nc is the number of unit cells in the perpendicular 
direction. This extra sum can be transformed into an 



J-Tr/2 TT ■ 



It is 



integral using the relation 
possible to perform the integrals analytically as we ex- 
plain in App. |B] Using the integrals defined there (/i 
and I2) we obtain for uj, Sa, Sb, t <C A 



rij 



-G 



AA 



A2 



-G 



Sa 



BB 



A2 



A2^ 



(37a) 
(37b) 



A2 



t±\uj\ 



e{2tj_ - \uj\ 



2A2 



1^1(1^1 + V^2_4i^) 



e{\uj\-2tjL), (39a) 



-^B - -^A - 



UJ 
A2 



sign(cij 



At 



A2 



ie(|w| - 2ti). 

(39b) 



The self-energies obtained within the t-matrix are shown 
in Fig. [T4l while those obtained from the CPA are shown 
in Fig. [m A comparison between the density of states 
in the different approximations is shown in Fig. 1161 In 
general the curves are similar to the ones in the bilayer 
but somewhat smoother. 

The behavior of the self-energies at the Dirac point in 
the multilayer are similar to the case of the bilayer treated 
in Section lV Al I V Bl and IV Cl We have more to say about 
this when we discuss the perpendicular transport in the 
multilayer in Section [IX Al 



B. Spectral function 

The spectral function for the graphene multilayer is 
given by a generalization of Eq. (|35p : 



^(k||,fcj^,t^) = -;^{lm[GgA(fc|hfcj^,^)] 

lm[Gg^{^,k^,uj)\]. (40) 



+ - 



Given the form of the Green's function and the CPA self- 
energies it is straightforward to obtain this quantity. The 
spectral function is depicted in Fig. [T7] for three values 
of the perpendicular momentum, since the model we use 
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FIG. 14: [color online] Self-energies within the t-matrix in the 
multilayer as a function of the frequency uj. Left: sublattice 
A; Right: sublattice B. 
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FIG. 15: [color online] Self-energies within the CPA in the 
multilayer as a function of the frequency uj. Left: sublattice 
A; Right: sublattice B. 



is electron-hole symmetric we only show the results for 
negative frequencies. We would like to stress that for a 
large part of the BZ the spectra are reminiscent of the 
bilayer spectra. At the edges of the BZ however, where 
2t± cos{k±d) = 0, the spectrum is that of massless Dirac 
fermions. 

Examples of the momentum distribution curves 
(MDC's) and the energy distribution curves (EDC's) in 
the disordered graphene multilayer are shown in Figs. fTSl 
and [19] for two different values of k± . The evolution of 
the peaks from delta functions to broader peaks with in- 
creasing disorder is clearly seen. One can also note that 
the influence of the impurities is more severe close to the 
H point of the BZ since the overlap of the peaks is larger 
there. The reason for this is that the particles there are 
dispersing linearly leading to peaks that are closer to- 
gether than for particles with a parabolic dispersion. 
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FIG. 16: [color online] Local density of states p on the differ- 
ent sublattices as a function of the frequency uj in the multi- 
layer. Top: sublattice A; Bottom: sublattice B. 
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FIG. 17: Intensity plot of the electron spectral function in 
the kxu plane for the multilayer for different disorder values 
and different values of k± in the CPA. From left to right: 
rii = 10"*, 10"^, 5 X 10"^. From top to bottom: fcx = 0, 
7r/(3d), 7r/(2d). 



VII. ELECTRON TRANSPORT IN BILAYER 
AND MULTILAYER GRAPHENE 



Having worked out the self-energies in the previous sec- 
tions we now turn to linear response (Kubo formula) to 
study electron transport. We saw in Sections IVl and IVTl 
that the low-energy states are mainly residing on the B 
sublattice. Nevertheless, electron transport coming from 
nearest neighbor hopping must go over the A sublattice. 
This is particularly important for the case of perpendicu- 
lar transport, since in the simple model that we are using, 
hopping comes exclusively from states with weight on the 
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FIG. 18; [color online] MDC's in the graphene multilayer for 
two values of the perpendicular momentum: fcx = (i.e. at 
the K point, parabolic dispersion, dashed line) and fcx = 7i"/2 
(i.e. at the H point, linear dispersion, solid line). The three 
panels are for different values of the density of impurities in 
the CPA. From the top the energy cuts are at the energies: 
.0005, .0055, .0105, .0155, .0205, .0255, and .0305 in units of 
the cutoff A. The curves are uniformly displaced for clarity. 
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FIG. 19: [color online] EDC's in the graphene multilayer for 
two values of the perpendicular momentum: k±_ = 7r/3 (i.e. 
between the K and H point, parabolic dispersion, dashed line) 
and k±_ = 7r/2 (i.e. at the H point, linear dispersion, solid 
line). The three panels are for different values of the density 
of impurities in the CPA. From the top the cuts are at fixed 
values of the in-plane momentum fc||: .0001, .01, .02, .03, .04, 
.05, and .06 in units of the cutoff A. The curves are uniformly 
displaced for clarity. 



A sublattice. 

This feature implies that although the total density 
of states at half-filling is finite, because the density of 
states on the A sublattice goes to zero as the Dirac point 
is approached, the in-plane and out-of-plane transport 
properties are unconventional. The main purpose of this 
section and the two following sections is to show how 
this comes about in detail through concrete calculations. 
We calculate conductivities (or optical response) parallel 



and perpendicular to the planes in both bilayer graphene 
and multilayer graphene. The resulting conductivities 
are very anisotropic and we find a universal nonzero mini- 
mum value for the in-plane DC conductivity as a function 
of the chemical potential. 



A. Conductivity formulas 

To calculate the conductivity we use the Kubo 
formulai^ We only consider the homogeneous (q ~ 0) 
response, but we investigate both the temperature de- 
pendence and the frequency dependence of the various 
conductivities. 

The conductivity is computed from the Kubo formula: 



cr(w) 



1 

Suj 



dte-\[j\tijm) 



uj -\- 15 



n(c.). 



(41) 

Here S is the area of the system, J is the current oper- 
ator of interest, and 11 the appropriate current-current 
correlation function. A contribution to 11 from a term of 
the form {[A{t), B{Q)\) where A = X)k '^('*^)'^ik'^2k and 
B = X^k /5('*)^ik^2k can be shown by the usual methods 
to give a contribution to the real part of the conductivity 
of the form: 



nF(e + cj) — n-p(t) 



k 

X Im[Gb2,ai(e,k)] Im[Ga2,w(e + k)] a(k)/3(k). (42) 

Here the imaginary parts only involve the frequency part 
and not the angular (spatial) parts of the propagators. 
In terms of the expression in Eq. ([^0]) this imply that the 
imaginary parts involves g'^ and g^'^ but not the spatial 
information encoded in Al]^(k) and Al|(k). With the 
inclusion of the two spin projections and two valleys we 
get (putting back h — 2'Kh and extracting the electric 
charge from the current operators): 



R.[.(„)] . ^ 



71^(6 + uj) — nF{e) 



(43) 

Here np is the Fermi distribution function. We have also 
introduced the kernel 5 that for the case of the operators 
above becomes: 



d0(k) 
2tt 



X Im[Gb2,ai(e,k)] Im[G,2,w(e + c^, k)] a(k)/3(k). (44) 

Thus the contribution to the in-plane DC conductivity 
at zero temperature is 

2e2 



nh 



S(0,0) = ao||S(0,0). 



(45) 



Finally we also note that - since we are using the approx- 
imation of purely local impurities - there are no vertex 
corrections appearing in the model. 
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B. Bilayer graphene 

The current operators can be obtained from the 
Peierls' substitution j^^i^^ and an expansion close to the 
K (or K') point in the BZ. Alternatively the current op- 
erators can be obtained directly from the Hamiltonian 
in Eq. ^ using J = env with the velocity being given 
by V = 9Ho(k)/9k. In any case, the current operators 
needed for the calculation of the conductivities are given 
by: 

Jxi = vpe [cAi,kCBi,k + cjsi.kCAi.k] > (46a) 

k 

Jx2 = Vpe [c^A2M^B2M + 42,kCA2,k] : (46b) 
k 

Jj_ = iedti_ Y [cAi,kCA2,k - c^2,kCAi,k] • (46c) 

k 

From the contributions of the form ([J^i, J^i]) to the 
current correlator we get a contribution to the kernel 
which is: 

^.i,.i= d(i;|fc2)|lm[gDA(e,k)]lm[ggB(e+^,k)] 
+ Im[ggB(e, k)] Im[gDA(e + ^, k)] }. (47) 
Similarly from the cross term ([Jxi, Jx2\) we get: 

.A- 

5.i,x2 = 2 / d{vlk^)lm[gll{e,]^)\lm[gll{e + u;,]^)\, 
Jo 

(48) 

while for the interplane optical response the contribution 
from ([J_L, J±]) is: 




X {lm[ggA(e,k)] Im[ggA(e + ^,k)] 

-Im[g^?(^^j^)]j^^gND(^_^^^j^)]y (49) 

Due to the phases in the Green's functions the other 
terms such as those involving G^gG^g vanish upon per- 
forming the angle average. To get the total CT|| per plane 
in the bilayer one should add the contributions from 
^xi,xi and ^xi,x2- 

C. Multilayer graphene 

The expressions for the current operators in the mul- 
tilayer are obtained in a similar way. Jy is given by the 
same expression as in Eq. (1461) except that the momentum 



variable is now three-dimensional. The current operator 
in the perpendicular direction is 

J± = -2et±d Y sin(A;_L)[c]^^ ,^CA2,k + c^2,kCAi,k]- 

k|| , k± 

(50) 

To get the conductivities in the multilayer, we should di- 
vide by the volume V = S x 2dNc instead of the area 
S. Here Nc is the number of unit cells in the perpen- 
dicular direction. We then turn the sums into integrals 
using ;^Efc^ IT- Thus to get cry we use the 

expressions in Eq. (|^ and Eq. (|47ll48p and add the per- 
pendicular integral according to: 




Similarly for the perpendicular conductivity we use 
Eq. gSl) with: 

X sin2(fc^){lm[gDA(e,k)] Im[g°A(e + ^, k)] 

+ Im[gSS(^^j^)]j^[gND(^_^^^j^)]|_ (52) 

The numerical value of the dimensionless prefactor in 
E!_L is approximately 0.15 (using d « 2.5a). When we 
present the results in the following sections it is con- 
venient to use a unit that combines the prefactor in 
this kernel with the factor from Eq. according to 
a^o = [2eyiTTh)]il/d)[Mt^/i3at)]\ 

VIII. RESULTS FOR THE CONDUCTIVITIES 
IN THE GRAPHENE BILAYER 

Using the formulas for the kernels (i.e., the S's) from 
the previous section and the explicit forms of the propa- 
gators from Section Hv] [in particular Eqs. ([20| -(|22 p ] we 
have calculated the kernels for arbitrary values of the 
self-energies. Details of this calculation are provided in 
App.O In this section we present results for the conduc- 
tivities [via Eq. (|43|) and Eq. (|45|) ] using the kernels ob- 
tained with the t-matrix and CPA self-energies discussed 
in Sec. El 



A. Chemical potential dependence 

The results for the DC conductivities at T = in the 
t-matrix and CPA approximations for different values of 
the chemical potential are shown in Figs. 1201 and [2T] The 
only difference between these figures are in the scales of 
the axes. It looks as if the minimum value for ay per plane 
in the bilayer is approximately 2cryQ = 4e^/(7r/i), which is 
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identical to the result one obtains using the same meth- 
ods in single layer graphene.^^ We discuss the minimum 
conductivity later in this section. 

The low-energy feature in the t-matrix curves comes 
about at the energy scale at which the two planes start 
to decouple. The scale at which this takes place (Sa ~ 
t±) is easily found numerically with the results shown 
in Table |TT1 The local maximum in the conductivities is 
readily identified with this energy scale. In the CPA this 



m 


0.01 


0.005 


0.001 


0.0001 


iu/A 


1.4 X 10"^ 


6.5 X 10"^ 


1.1 X 10"^ 


8.4 X lO"-"^ 



TABLE 11: Energy scale cu at which the planes start to de- 
couple within the t-matrix approximation [this happens when 
T.a{co) ~ tj_]. 

scale is suppressed and the curves show no peak. Quite 
generally it is seen that the CPA curves are smoother 
than the ones for the t-matrix. 

Another interesting feature of a_L is that it is increased 
by disorder. This is due to the fact that disorder enhance 
the DOS on sublattice A where all the contribution to 
(Tj_ is coming from. At the lowest values of the chemical 
potential the perpendicular conductivity still goes to zero 
however. 



B. Minimal DC conductivity 

By studying the curves more closely, it looks as if the 
CPA curves actually gives a value that is smaller than 2 
in the limit lu 0. In fact, we find that the minimum 
in the in-plane DC conductivity is again (as in the sin- 
gle layer case of massless Dirac fermions^-) universal in 
the sense of being independent of the particular impu- 
rity concentration. In the bilayer the minimum value per 
plane obtained from the CPA is 



f||min,CPA — " T" 



3e^ 

TT h 



(53) 



This value is obtained by using the form of the self- 
energies in the low frequency limit that are given in 
Eq. (P^ . Explicitly one finds for the propagators via 
Eq. jUl): 



gAA(w 



gBB('^ 



CA2 + fc2 ■ 
CA2 + fc2 ' 

0. 



(54a) 

(54b) 
(54c) 



>0,k) 

>0,k) 

g^g(c.->0,k) 

Using these asymptotic forms in Eq. (|47|) and Eq. 
the contribution from the latter equation drops out. The 
value in Eq. (|53p is obtained from the first term after em- 
ploying the relation: ImpAlImpe] (\/3/2)2|I;aEb| 
3eA2/4. 



We note that our value for the minimal conductivity 
is different from the values obtained in works by other 
groups. In particular both Koshino and Ando (using a 
2-band model in conjunction with a second-order self- 
consistent Boltzmann approximation) and Snyman and 
Beenakker (using the conductance formula for coherent 
transport) both finds the value Ae'^/{TTh) per plane j^2i57 
(The minimal conductivity problem in bilayer graphene 
has also been discussed in Ref. [ssllsol ). We can repro- 
duce their result in our formalism by considering the case 
that the impurities only affect the B sublattice, as dis- 
cussed in Sec. IV Bl (or the case of Born scattering dis- 
cussed in Sec. IV C| . In particular, taking Sa = and 
Sb(w = 0) = -iT [from Eq. dST])] one finds 



?gA(^-0,k) 
?Ag(^ = 0,k) 



(Ur)2 + A:4' 
-iTtl 

(iir)2 + ' 

iTt±k 



(55a) 
(55b) 
(55c) 



Using these expressions in Eq. (|47p and Eq. i|^o|) . ^xi,xi 
and '^xi.x2 are found to contribute equally to the con- 
ductivity with the total value being 4e2/(7r/i). This re- 
sult shows that the minimal conductivity is not really 
"universal" in the graphene bilayer since it actually de- 
pends on how the impurities are distributed among the 
inequivalent sites of the problem. Furthermore, the bal- 
listic results corresponds to the case where the disorder is 
only affecting the B sublattice. Nevertheless, the general 
conclusion is that there is a non-zero minimum in the in- 
plane conductivity of the order of e^//i. Further evidence 
for this conclusion is hinted by the related issue of how 
other hopping terms, such as 73 and 74, affect the value 
of the minimal conductivity. The case of trigonal warp- 
ing (i.e. 73 ^ 0) have been discussed in detail by Cserti 
and coUaborators,^^^ and they find that the minimal value 
of the conductivity per plane is {12 / tt) {e'^ / h) in this case 
(See App. I^for an alternative derivation of this result). 
The introduction of 74 (or A, or a next-nearest neighbor 
hopping within the planes) - which breaks the symme- 
try in energy between the central Dirac point and the 
three elliptical cones away from k = will likely further 
influence the minimal value. 



C. Frequency dependence 

The frequency dependence of the conductivities are 
pictured in Fig. [22ll25l The figures reveals some inter- 
esting features of the band structure and the semimetal- 
lic behavior. For the case of a finite chemical potential 
the temperature does not make a big difference since it 
(at 300 K) is still much smaller than the Fermi energy, 
this is manifested in the small difference between Fig. [Ml 
and Fig. [551 Near zero chemical potential the effect of 



16 



T-matrix, T = 0, w = 




A 



T-matrix, T = 0, = 



n, = 0.0001 

n, = 0.001 

n, = 0.005 



!/>=-.. 



0.05 0.1 0.15 0.2 0.05 0.1 0.15 0.2 




FIG. 20: [color online] DC conductivities in the bilayer as a 
function of the chemical potential (in units of the cutoff) at 
zero temperature. Left: t-matrix; Right: CPA. Top: in plane; 
Bottom: c-axis. 



FIG. 22: [color online] Conductivity as a function of frequency 
(in units of the cut-off) for the bilayer at the Dirac point 
= for T = 0. Left: t-matrix; Right: CPA. Top: in plane; 
Bottom: c-axis. 
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FIG. 21: [color online] DC conductivities in the bilayer as a 
function of the chemical potential (in units of the cutoff) at 
zero temperature. Left: t-matrix; Right: CPA. Top: in plane; 
Bottom: c-axis. 



FIG. 23: [color online] Conductivity as a function of frequency 
(in units of the cut-off) for the bilayer at finite temperature 
T = 300-?sr near the Dirac point = 0. Left: t-matrix; Right: 
CPA. Top: in plane; Bottom: c-axis. 



the temperature is more dramatic. The temperature in- 
crease the number of carriers and is responsible for the 
Drude-like peaks that appear in the plots for low impu- 
rity concentrations. A well-known feature of semimetals 
is that the temperature is an important factor in deter- 
mining the number of carriers in the system. The peak 
at w ^ = .05A is due to the onset of interband 
transitions. The frequency-dependence of the absorp- 
tion of electromagnetic radiation has also been studied 
by Abergel and Falko with similar results^^ they also 
study transitions between Landau levels in a magnetic 
field. 



D. Temperature dependence 

The temperature dependence of the DC conductivity 
can be found in Fig. and Fig.[271 For the case of a finite 



chemical potential the in-plane conductivity curves are 
flat and proportional to l/n^, as is expected in a normal 
Fermi liquid metali^ The contribution to the scattering 
from impurities is very weakly temperature dependent. 
Nevertheless, there is a small temperature dependence for 
the lowest impurity concentration which is due to the fact 
that T/E-p is still not negligible. Near the Dirac point the 
characteristics of a semimetal appear again as the con- 
ductivities become temperature dependent. Note how- 
ever that we are not considering scattering by phonons 
which is important at finite temperatures. 



IX. RESULTS FOR THE CONDUCTIVITIES IN 
MULTILAYER GRAPHENE 

Using the same procedure as for the bilayer in the pre- 
vious section we have calculated the kernels for arbitrary 
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FIG. 24: [color online] Conductivity as a function of frequency 
(in units of the cut-off) for the bilayer at finite chemical po- 
tential fi = 0.025A at T = 0. Left: t-matrix; Right: CPA. 
Top: in plane; Bottom: c-axis. 



FIG. 26: [color online] Temperature dependence of the DC 
conductivities in the bilayer at the Dirac point fi — 0. Left: 
t-matrix; Right: CPA. Top: in plane; Bottom: c-axis. 
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FIG. 25: [color online] Conductivity as a function of frequency 
(in units of the cut-off) for the bilayer at finite temperature 
T — 300K at finite chemical potential fi — 0.025A. Left: 
t-matrix; Right: CPA. Top: in plane; Bottom: c-axis. 



FIG. 27: [color online] Temperature dependence of the DC 
conductivities in the bilayer at finite chemical potential /i — 
0.025A. 



values of the self-energies. Details of this rather lengthy 
calculation are provided in App. [Cl In this section we 
present results for the conductivities using the kernels 
obtained with t-matrix and CPA self-energies discussed 
in Sec. ED 

The DC conductivities in the multilayer as a function 
of the chemical potential fi are pictured in Figs. [^511^ 
the only difference between the figures are in the scales 
of the axes. The property of disorder-enhanced transport 
in the perpendicular direction seems to survive in this 
model for the multilayer, but only for very low values of 
the chemical potential. For larger values of the chemical 
potential the influence of disorder becomes more conven- 
tional. In this case, because of the finite Fermi surface, 
the transport properties are more like in a normal metal. 



A. Perpendicular transport near the Dirac point 

Generalizing Section IV CI to the multilayer we find 
again that Sa ~ and Sb ^ — iFe is purely imaginary 
in the Born limit at the Dirac point. Nevertheless, as we 
shall see it is necessary for the computation of (j± that 
Sa remain finite. Therefore we take Ea ~ —i^A and as- 
sume that Fa ^ Fb . We note that this is also consistent 
with a self-consistent version of the Born approximation 
for weak potentials. Thus, for w — s- we have 



T [ D 1 -FB(fc^+FAFB) . 

^"^^■^^^ ^ [2<^FBCOs(fc^)]2 + (fc2+rAFB)2^^^^^ 

MsIa] - 0. (56b) 
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Inserting these expressions into Eq. ((52)) it is possible to 
perform the integrals exactly with the result 



^ J_,inulti 



1 Mdt± 
d V Sat 



log 



^1 + (rA/2ti)2 + 1) 



Vi + (rA/2ti)2-i) 



(57) 



Thus there is a logarithmic singularity in the limit Fa — > 
as mentioned above. Intuitively this singularity comes 
from "clean" chains of atoms along the A sublattice 
where transport is unhindered once some weight has been 
pushed onto the A sublattice by the impurities on the B 
sublattice. It is plausible that S_L_niuiti increases with in- 
creasing disorder. It is so because the first factor grows 
linearly whereas the second factor decays only logarith- 
mically with the r in question. 

For the case of vacancies in the CPA a result analo- 
gous to the one in Eq. ([29]) can be obtained. In fact the 
result is the same up to a factor: Ea — > 2^/^Sa and 
Eb 2~"'^/'^Eb- Therefore the asymptotic expressions 
in Eq. (j54p are valid also in the multilayer. In addition 
one finds that gAAl"^ ^ 0, k) 0. Thus, asymptotically 
one finds that S±,muiti ^ w^/"^, which leads to a temper- 
ature dependence of a± at the Dirac point that is of the 
form T^/^. We also note that 

f||,min is independent of 
t± in the bilayer, thus we conclude that it takes on the 
same value in both the bilayer and the multilayer. Using 



the fact that S 



II ,niulti 



constant at the Dirac point we 



find that a\\/aj_ diverges as T as T 
previouslyiii 



as reported 
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FIG. 28: [color online] DC conductivities in the multilayer as 
a function of tiie chemical potential (in units of the cutoff) 
at T=0. Left: t-matrix; Right: CPA. Top: in plane; Bottom: 
c-axis. 



B. Frequency and temperature dependence 

The frequency dependence of the conductivities in the 
multilayer is shown in Figs. [5011331 for two different tem- 
peratures and both at the Dirac point and for a finite 
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FIG. 29: [color online] DC conductivities in the multilayer as 
a function of the chemical potential (in units of the cutoff) 
at T=0. Left: t-matrix; Right: CPA. Top: in plane; Bottom: 
c-axis. 
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FIG. 30: [color online] Conductivity as a function of frequency 
(in units of the cutoff) at the Dirac point /x = for T = in 
the multilayer. Left: t-matrix; Right: CPA. Top: in plane; 
Bottom: c-ajds. 



chemical potential. For the cleaner systems a Drude-like 
peak appears at finite temperatures for both in-plane and 
perpendicular transport at the Dirac point. For a finite 
chemical potential - because the system is metallic in 
both directions - the system has a Drude peak in the 
conductivity also at zero temperature. Moreover, it can 
be seen how the suppression of the conductivity in the fre- 
quency range before interband contributions sets in (i.e. 
at w = 2/i) is affected by both disorder and temperature. 

We note that our curves for the frequency-dependent 
in-plane conductivity is very similar to the recent results 
of Ref . "e^, which include both measurements and calcu- 
lations based on the full SWM model. 

Our results for the temperature dependence of the con- 
ductivities in the multilayer are shown in Figs. [34ll35l At 
the Dirac point, the in-plane conductivity goes to a fi- 
nite constant while the perpendicular conductivity goes 
to zero as T — > 0. The disorder-enhanced transport at 
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T-matrix. T = 300K, n = 0.025 A 
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FIG. 31: [color online] Conductivity as a function of frequency 
(in units of the cutoff) at finite temperature T — 300K near 
the Dirac point fj, — in the multilayer. Left: t-matrix; Right: 
CPA. Top: in plane; Bottom: c-axis. 



FIG. 33: [color online] Conductivity as a function of frequency 
(in units of the cutoff) at finite temperature T = 3007^ at 
finite chemical potential fi = 0.025A in the multilayer. Left: 
t-matrix; Right: CPA. Top: in plane; Bottom: c-axis. 
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FIG. 32: [color online] Conductivity as a function of frequency 
(in units of the cutoff) at finite chemical potential jj, — 0.025A 
at r = in the multilayer. Left: t-matrix; Right: CPA. Top: 
in plane; Bottom: c-axis. 
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FIG. 34: [color online] Temperature dependence of the DC 
conductivities at the Dirac point fi = in the multilayer. 
Left: t-matrix; Right: CPA. Top: in plane; Bottom: c-axis. 



low temperatures can also be seen in the figure. For a 
finite chemical potential, the system behaves like a metal 
with only a weak temperature dependence of the conduc- 
tivities. 



present here were previously reported in a brief form in 
Ref. [l3|- A recent detailed study of the impurity states 
in the unitary limit in both biased and unbiased bilayer 
graphene can be found in Ref. [ilj. 



X. IMPURITIES IN THE BIASED GRAPHENE 
BILAYER 



In the following sections we study the problem of im- 
purities and mid-gap states in a biased graphene bilayer. 
We show that the properties of the bound states, such 
as localization lengths and binding energies, can be con- 
trolled externally by an electric field effect. Moreover, 
the band gap is renormalized and impurity bands are cre- 
ated at finite impurity concentrations. Using the CPA we 
calculate the electronic density of states and its depen- 
dence on the applied bias voltage. Many of the results we 



A. Band model 

In this section, we review the properties of the minimal 
model introduced in Eq. PH]). Throughout this section 
we use units such that vp = Ti = 1 unless otherwise 
specified. For numerical estimates we use t± = .35 eV 
and insert the appropriate factors of v-p — Ma/2 with 
t = 3eV and a = 1.42 A. From Eq. (UHl) one finds two 
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FIG. 35: [color online] Temperature dependence of the DC 
conductivities at finite chemical potential fi — 0.025A in the 
multilayer. Left: t-matrix; Right: CPA. Top: in plane; Bot- 
tom: c-axis. 



pairs of electron-hole symmetric eigenvalues: 



4 

_L' 



(58) 



where s — zk. The lowest energy bands (with respect to 
the "Dirac point" at zero energy) representing the valence 
and conduction bands are the -E±,- bands. The smallest 
gap takes place at a finite wave vector given by 



V I V^ + 2tl 
" " 2 V V^+tl' 

so that the size of the band gap becomes 



(59) 



At k = the distance between the valence and the 
conduction band is given by the applied voltage differ- 
ence V. Note that V should in reality be taken to be 
not the bare applied voltage difference but instead the 
self-consistently determined value Vmf as discussed in 
Refs. [llHTslflsHs^ l . Near fcg the energy of the quasi- 
particles in the conduction band can be expanded as 



E+.- 



2tl) 



4^3/2 



2m* 



(61) 



and as long as this approximation is valid the density of 
states per unit area is 



N{uj) 



k() 



2m* 



(62) 



^ Y l^l-Sg/2' 
for \uj\ > Eg/2. This includes both the valley and the spin 
degeneracy. Notice that the divergence of the density of 
states (DOS) at the band edge is similar to what one 
would get in a truly ID system. The fact that a large 
DOS is accumulated near the band edge has important 
consequences for the properties of the bound states as we 
shall see in the following. 



B. Bare Green's function 



An explicit expression for the bare Green's function, 

1 . 



(60) which is given hy = [uj ~ Hbb] , is: 



D 



f {uj- V/2) [{uj + y/2)2 - fc2] [(^cj + Vl2f - fc2] ke^'i' 

{{lo + y/2)2 - fc2] _ -^/2) [(w + Vl2f - e\ ~{uj + V/2)t\ 

t^{Lj^-V^/A) t^{Lu + V/2)ke"^ 

\tj_{u;~V/2)ke"l' t^k^e^'"!' 

tj_ (cj2 - V^/i) tj_ (w - V/2)ke-'"t' 

t^{uj + Vl2)ke-^'t' t_Lfc^e-2# 

{uj + V/2) [{uj - V/2f - e] [{lu - V/2f - e] ke-"*' 

[{uj - V/2f - fce*-^ {uj + V/2) {{uj - V/2Y - fc^] 



, (63) 



V 



where 



D 



VH^ 



■uj\V^^t\). 



(64) 
(65) 



So that, for example, the important diagonal components 
are given by 



G 



AlAl 



{uj-V/2){{io + V/2f -k'' 



po _ po {^^Y/2)t\ 



D 



(66a) 
(66b) 
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The corresponding components for plane 2 are obtained 
by the substitution V —V. Note that M > inside of 
the gap. 



XI. BOUND STATES FOR DIRAC DELTA 
POTENTIALS 



As mentioned previously, to locate the bound states we 
must find possible new poles due to the potential. Ex- 
plicitly we need UG (e) = 1 . Like in Section IVl G is 
the local propagator at the impurity site that is given by 
the expression in Eq. ([25ll with G° taken from Eq. (|63|) . 
Using Eq. ()66p we can perform the integrals exactly in 
the continuum approximation with the result 



Bound states must be located inside of the gap so that 
their energies fulfill |e| < Eg/2, otherwise the asymptotic 
states at infinity are not exponentially localized. If we 
decode a number (say Ni) of local impurities in a matrix 
of the form 



\> = Diag[f/i, U2,...,Un,], 



(67) 



where we let Ui denote the strength of the impurity po- 
tential that is located at site x^. The total Green's func- 
tion is then given by 



G 



G° + gVg 
= G° + G"[^/- 



: G° + gVg° + G°VG"VG° + ... 
VG^V + VG^VG'^V + . . .]G° 

EE G° + G°TG°. (68) 



Here T is the t-matrix of the system (see e.g. Ref. [5l|V 
The interpretation of this expression is most transparent 
in the real space picture, where it includes the repeated 
scattering off of all of the impurities in every possible 
way. Another way of expressing T is (decomposing V as 



T= VV{1 



VG" 



V + i\/VGWVf + ...)\/v 



1 



fgVf 



(69) 



Bound states generated by the impurities can readily be 
identified by the possibly new poles in the full propagator 
of the system. Therefore an equation that can be solved 
to find the energies of the bound states of the system is 
given by 



Det 



C/,;G°-(6) 



0. 



(70) 



Here G°j (e) denotes the (real space) propagator from site 
J to site i. In principle one can put in an arbitrary num- 
ber of impurities in this expression. However, if two im- 
purities are located too close to each other the continuum 
approximation to the propagators is not expected to be 
accurate and one must instead work with the full tight- 
binding description (see Ref. [65| for an illustration of 
this approach in monolayer graphene). If we specialize 
to one local impurity affecting only one site the calcula- 
tions simplify considerably. The Fourier transform of the 
local potential \sU/N (where N is the number of unit 
cells in the system) so that we can write 



T 



U 



N 



C/G°' 



(71) 
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(72b) 



where A'P = ^VH\/ A - uj'^{V'^ + t\), and A (w 7 eV) is 
the high energy cutoff. The corresponding expressions in 
plane 2 are obtained by the substitution V ^ —V. The 

typical behavior of G (co) as a function of the frequency 
uj is shown in Fig. 1361 
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FIG. 36: [color online] Left: typical behavior of the G 's: 

Gbi{(^) (dotted), Gb2('^) (dashed), Gai(w) (solid), Ga2('^) 
(dash-dotted). Here V — 50meV. The divergences close to 
the band edges are clearly visible. Right: bound state energies 
e for a local potential of strength —U, the labeling of the 
sublattices is the same as to the left. 

From this we conclude that a Dirac delta potential al- 
ways generates a bound state (no matter how weak the 

potential is) since G diverges as the band edge is ap- 
proached (where ill — > 0). The dependence on the cutoff 
(except for the overall scale) is rather weak so that the 
linear in-plane approximation to the spectrum should be 
a good approximation as in the case of graphenci^ For a 
given strength of the potential U, there are four different 
bound state energies depending on which lattice site it 
is sitting on. In Fig. [321 we show the energies of these 
bound states for strong impurities. Even at these scales 

the bound state coming from G^i is so weakly bound that 
it is barely visible in the figure. In Fig. [37] we show the 
binding energy as a function of U and V for the deepest 

bound state (coming from Ggj^). In the limit of — > cxd 
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the electron-hole symmetry of the bound state energies 
is restored. For illustrative purposes we consider only 



V = 200meV 

0.025 V= 100 meV 

V = 40 meV 
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FIG. 37: [color online] Left: Bound state binding energies 
iSb (in units of the gap E^) for a Dirac delta potential of 
strength —U for different bias V. Right: Binding energies of 
a potential well of range R = 10 a and strength 7 = 71 = 72 
(see the text in Sec. IXIV)| for different angular momentum m 
channels and external bias potential V. 

attractive potentials in this work, analogous results hold 
for repulsive potentials because of the electron-hole sym- 
metry of the model that we are using. For smaller values 
of the potential {\U\ ^ A) the binding energy measured 
from the band edge E-^, — Eg/2 — e grows as C/^ and the 
states are weakly bound. For example, for V = 40meV 
and U < leV one finds < 4 x IQ-^E^. 



A. Angular momentum states 

For any potential with cylindrical symmetry it is use- 
ful to classify the eigenstates according to their angular 
momentum m. In the presence of the "trigonal distor- 
tion" parametrized by 73 the calculations become more 
involved because of the broken cylindrical symmetry. We 
discuss this issue briefly in Sec. lXVll The real space con- 
tinuum version of the Hamiltonian matrix in Eq. (I10|) 
that includes a potential, which in general is allowed to 



take on different values in the two planes, is 



Ho = 



_ V/2 + gi{r) -ia* -d 



t^il + (7,)/2 

-V/2 + g2{r) - ia ■ d . 



(73) 

Here ai {i — x,y,z) are the usual Pauli matrices. For 
example, a symmetric Coulomb problem could then be 
approximated by taking gi{r) = g2{j') = gjr. Going to 
cylindrical coordinates the derivatives transforms accord- 
ing to 



dx ± idy — e 



(74) 



where we use the coordinate convention x ± iy 



For the Hamiltonian in Eq. (|73p one can now - in analogy 
with the usual Dirac equation^® - construct an angular 
momentum operator that commutes with the Hamilto- 
nian. The angular (ip) dependence of the angular mo- 
mentum m eigenstates are those of the vector: 



/ 1 \ 

1 



V 



(75) 



where parameter a is introduced for later convenience, it 
is used later to obtain more compact expressions. It is 
convenient to define the following "star" product of two 
vectors that results in another vector with components 
given by 



(76) 



By writing 4* = uo,m * tpif)/ the eigenvalue problem 
Hq"^ = E'i> is transformed into a set of coupled ordinary 
linear differential equations for the radial wave-function 

ijj{r): 



( gi{r) + V/2 -idr + ij/r 

-idr - ij/r gi{r) + V/2 
t± 
\ 



t± 



92{r) - V/2 
-idr + i (j + l)/?" 



\ 


'idr -i{j + l)/r 
.92(0 - V/2 ) 



ijj{r) = Eip{r) 



(77) 



Here we have introduced j = m — 1/2 to render the equa- 
tions more symmetric. If the potential generates bound 
states inside of the gap these states decay exponentially 
~ r'^'e^'"' as r — > 00. Assuming that the potential decays 
fast enough the asymptotic behavior of Eq. (|77p imply 



that the allowed values for k are k± satisfying 

K± = ^/-{e^ + V^/A)±iM^, (78a) 
|k|4 ^ (FV4-e')(FV4-e'+ii), (78b) 

K± = |At|exp{±4|-itan-i(-^^^;^)]}(.78c) 
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So that, for weakly bound states we have 



C. Local impurity wave functions 



.1 



K± S 

leading to a localization length 

^kg- I Eg; 



±i-^V^ + El, (79) 
1/2 + K2 2 




(80) 



that diverges as the band edge is approached and de- 
creases as the bias voltage increases. 



B. Free particle wave functions in the angular 
momentum basis 



The general expression for the retarded Green's func- 
tion is 

OJ + 10 — hn 

n 

where the sum is over the eigenstates \n) (with eigenen- 
ergy En) of the system. Comparing the coefficient of the 
poles in this expression with those in Eq. one can 
read of the wave functions of the bound states directly. 
The result is that the wave functions are Fourier trans- 
forms of the columns of the bare Green's evaluated with 
the frequency set to be equal to the energy of the bound 
state. 



The free particle wave functions in the angular mo- 
mentum basis can be conveniently expressed in terms of 
the following vectors: 



'fZ,m{z) 



' Zrr,{z) \ 

Zm{z) 
\Z„i+i(z)/ 



(81) 



w{p) 



( [(w + V/2f - p2] (cu - V/2)\ 
tj_{uj - V/2)p 



(82) 



The last vector is useful as long as ^ V/2 (cf. the 
discussion of the two eigenvectors in Ref. [11].) The de- 
nominator (actually a determinant) that determines the 
eigenstates is: 



D{p,L0) 



P' 



2 VyA-u;^f+Vhl/4-u;\V^+tl). (83) 



Then, provided that D{k,uj) = 0, (fc > 0) which cor- 
responds to propagating modes, the eigenfunctions are 
proportional to: 

*Z,m(t^,fc,f) = Ul,„i *WZ,m(^7') *w(fc), (84) 

where Zm{x) = Jm{x) or Ym{x) are Bessel functions and 
the star product is defined in Eq. (ffS)) . If on the other 
hand D{iK,uj) — 0, (Re[K] > 0) the eigenfunctions are 
instead: 

*K,m(w,K,r) = UQ,m'*^VK,TniKr)'kw{iK), (85) 

5'/,m(w,K,r) = wq,™ ★t>/,™(Kr) ★ w(-iK), (86) 

with Im{x) and Km{x) being modified Bessel functions. 
That these vectors are indeed free-particle eigenstates 
can be verified straightforwardly by applying the Hamil- 
tonian in Eq. (|73p with 51 = (?2 = to them. 



1. Impurity on an Al site 

When the impurity is on an Al site the expression 
becomes [using Eq. 
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Performing the angular average one ends up with Bessel 
functions: 



f Gaiai\ 
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\ it^{uj -V/2)kJi{kr)e^'f 

(89) 

There are really two such terms, one for each K-point, 
which corresponds to the different signs of the phases 
in Eq. ([55)1 . Note that this state has angular momentum 
m = in the language of the Section lXI Al The fc- integral 
can be performed analytically (taking A ^ 00 in the 
integration limit) . Using k,± defined in Eq. ([78]) we obtain 
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These propagators can also be easily expressed in terms 
of the free-particle wave functions: 
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(91) 



24 



This property is not a coincidence since the particles are 
essentially free, except for the potential that acts on the 
single impurity site at the origin. 



2. Impurity on a Bl site 

When there is an impurity on a Bl site one can per- 
form the same calculation with the result that the wave 
function becomes 



G 



ai,Bl 



2A/2A2(F/2 



(92) 

This state has angular momentum m = 1 in the language 
of Sec. IXI Al A similar expression can be obtained when 
there is an impurity on an A2 (B2) site, where in this case 
the corresponding state has angular momentum m = 
(m = -1). 



D. Asymptotic behavior 

The asymptotic behavior of the modified Bessel func- 
tions is Kn{z) ~ exp(— 2)/-y/z as z — > oo. Therefore the 
bound states are exponentially localized within a length 
scale given by [See Eq. ([75]) ] 
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(93) 



where the last line is applicable for weakly bound states 
close to the band edge. This is in agreement with the 
general results above in Eq. ([50]) . At short distances one 
may use that Kn{z) ~ l/z" for n > and Ko{z) ^ — Inz 
to conclude that the wave functions are not normalizable 
in the continuum. In particular, for an impurity on the 
Al (Bl) site the wave function on the Bl (Al) site di- 
verges as 1/r. This divergence is however rather weak 
(i.e., logarithmic) and not real since in a proper treat- 
ment of the short-distance physics, the divergence is cut- 
off by the lattice spacing a (this is equivalent to cutting 
of the fc-integral in Eq. (|89p at fc = A instead of taking 
A oo). 



XII. SIMPLE CRITERION 

Using the asymptotic form of the wave functions one 
can approximate the wave function as: 



in each plane. Normalization then requires that A — 
k'/tt, where k! is the real part of k±. Thus one impu- 
rity is interacting with approximately 



Ni = (Trr^)/ 



3V3a2 



(95) 



atoms per plane. For an impurity density of rXj, the 
number of impurities interacting with a given impurity 
is given by 



(96) 



A simple estimate of the critical density ric above which 
the interaction between different impurities are impor- 
tant is then given by iV^ '--^ 1. Writing e ~ Eg/2 — E^, one 
then finds the following criterion for overlap of impurity 
wave functions (assuming weakly binding impurities): 



ni 



> 



1 



27rV3 ^ kgt ^ Eg 



2.5 X 10" 



(97) 



indicating that the critical density increases with the ap- 
plied gate voltage. The last step is valid for V t±. 
Taking U < leV we found in Section |Xl] that Eh < 



4 X 10^'^Eg, leading to Uc ~ 10~^. Hence, even tiny con- 
centrations of impurities lead to wave function overlap. 
This result shows that even a small amount of impurities 
can have strong effects in the electronic properties of the 
BGB. 



XIII. VARIATIONAL CALCULATIONS 

For general potentials it is not possible to solve for the 
bound states in closed form. Nevertheless, for estimates 
and to gain intuition about the problem it is fruitful to 
study the problem with variational techniques. In this 
section we consider two different variational approaches. 

A. Variational calculation I 

Using Eq. (1771) can show the existence of bound 
states variationally. For simplicity we consider only the 
case m = {j ~ —1/2) and a symmetric potential gi — 
92 — gir)- We use the simple trial wave function ip2 — 
v4exp(— A:r/2). The following integrals are useful in the 
process: 



1, 



(98a) 





oo 



1^2 



-dr 



Q{R-r)\il:2\^dr 



kEi{ka), (98b) 
1 - exp(-fci?). (98c) 



(94) 



Here a is a cutoff on the order of the lattice spac- 
ing needed to regularize the integral. Ei{x) = 
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exp(— r)/r is an exponential integral (see e.g. 
Ref. [l^l). The eigenvalues (ck) of the kinetic term is 
given by the equation 



= {e^ + ^-fL[E,ikaf^l]} 



(99) 

which is the same as the equation for the bare bands 
[cf. Eq. jMl)] with the substitution k"^ k'^[Ei{ka)'^ - 
l]/4. Provided that ka < 0.26 (fc < 1.2 eV) the new 
"momentum" is real. For an attractive potential we may 
then construct a wave packet corresponding to the 
band leading to a positive contribution from the kinetic 
term. 

We consider two types of potentials: one of the 
Coulomb type, gc — —otjr, characterized by the di- 
mensionless strength a; and a local potential, gi^ = 
—{/(^{R — r), characterized by the strength U and the 
range R. The total variational energies for the two types 
of potentials are: 



Evar,L = E+- 



VEiiko 



1 



akEi{ka)(,100a.) 





FIG. 39: [color online] Variational energy for a short-range 
potential of range _R as a function of the variational parameter 
k ioT V = 50meV. From top to bottom: U = .033, .1, .33, 
and 1 eV. Left: R=10a; Right: R=la. 



B. Variational calculation II 

Another simple variational approach is to construct a 
wave packet with angular momentum m and a momen- 
tum close to kg from the band using the free particle 

solutions in Eq. ([84]) according to: 



y/Ei{kaf - l) -U[l- e-^fpOb) *var(0 



dp'^j,,n{E+_ip),p,r) 



p 

4<' 



(101) 



Some typical results obtained from these expression are 
shown in Figs. [551 and [551 For a sufficiently strong poten- 
tial it is favorable for the state to become very localized 
close to the impurity, and the assumed "bound state" 
is located inside of the continuum of the valence band. 
This is problematic as it leads to a breakdown of the 
picture of a bound state coming only from the states in 
the band. The state can no longer be considered 

a true bound state since it is allowed to hybridize with 
the states in the valence band and hence leak away into 
infinity. This state can therefore only be regarded as a 
resonance. Nevertheless, for weak Coulomb potentials 
there are indeed bound states inside of the gap, and for 
short-range potentials the variational treatment give re- 
sults that are consistent with the more rigorous study 
coming up in Section [XIVI 




FIG. 38: [color online] Variational energy for a Coulomb 
potential as a function of the variational parameter k for 
V = 50meV. From top to bottom: a — .033, .1, .33, and 
1. Left: large view; Right: zoom in for small k. 



where we assume that ^ <C fcg. The factor \/p/£, is 
included to generate a properly normalized variational 
state. Since this state is built up of eigenstates of the 
kinetic term the contribution from the kinetic energy to 
the variational energy becomes: 



E. 



■.kin ~ >. 



dp 
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2m* 
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(102) 



using Eq. ((6T|) . The leading contribution to the interac- 
tion energy for small ^ becomes: 



var.iiit 



dp I dp 



X ^^l,jE+^{p),p,r)gir)^j,r,.iE+^ip'),p',r) 



2ik^ j drr-^\^{E, fcg, r)g{r)-^ j^m{E, /cg, r) 



iJ=B+_(feg) 

^-2C/cg-^. (103) 

Therefore, the variational calculation shows that for any 
m, a weak attractive potential of strength oc U leads to a 
weakly bound state with binding energy E^, oc . This 
can be understood by noting that for each value of m, 
due to fcg being nonzero, the problem maps into a ID 
system with an effective local potential. It is well known 
(see e.g. Ref. Q) that in ID a weak attractive potential 
(c!C U) always leads to a bound state with binding energy 
E\j (X U'^. Thus the result is a direct consequence of the 
peculiar topology of the BGB band edge - see however 
Sec. [XVll 



26 



XIV. POTENTIAL WELL 

For the case of the simple local "potential well" defined 
by the potentials — ~Uj<d{R~r) = —jjQ{R—r)/Rit 
is possible to make analytic progress with the continuum 
problem. In Sec. IXI Bl we gave the explicit form of the 
eigenstates for a constant potential in the angular mo- 
mentum basis. Bound states are possible when the two 
solutions for r < R and the two solutions for r > R are 
not linearly independent at r = i?. This can be tested by 
evaluating the 4x4 determinant of the matrix built up 
by the four eigenstate spinors. Given Ui, U2 and R the 

I 

Do{lu) 

D2{UJ) 



depending on whether there are zero, one or two propa- 
gating modes at the chosen energy inside of the potential 
region. Here, 

K± = V-(^^ + ^V4)±iM2, (105) 
p± = \j {Zj-^ + V'^/4)±\I -M^, (106) 

= \j\l-M'^-{u)^ + V'^/A), (107) 

where M is given by Eq. ([55)1 with the substitutions V 
V and uj ^uj. 

By monitoring the zeros of -D„ as a function of the ra- 
dius R and the strengths 7^ we have studied the binding 
energies and find that the deepest bound states are in 
one of the angular momentum channels to = 0,±1 for a 
substantial parameter range. Since these types of states 
are also present for the Dirac delta potential we argue 
that the physics of short-range potentials can be approx- 
imated (except for the short-distance physics) by Dirac 
delta potentials with a strength tuned to give the correct 
binding energy. A typical result for the binding energies 
is shown in Fig. [371 

A feature of potentials with a finite range is that upon 
increasing the potential strength the binding energies can 
be made to increase until the state merges with the con- 
tinuum of the lower band and becomes a resonance. This 
is illustrated in Fig. l40l where we have plotted D2{oj) for 
different values of the strength of the potential; and it 
can be seen how the zeros of D2 moves across the gap 
and ultimately disappears into the valence band. Notice 
that this is consistent with the interpretation of the vari- 
ational calculation of Sec. IXIII Al We expect a similar 
behavior to occur for a strong Coulomb potential, but 
this interesting case is beyond the scope of this study. 



resulting determinant is a function of the energy w. Ze- 
ros of the determinant inside of the band gap corresponds 
to the bound states that we are searching for. Inside of 
the potential region the effective frequency and bias are 
given by: 

S = CO + (t/i + C/2)/2, (104a) 
V = V + {U2- Ui). (104b) 

The determinant is given by one of the following expres- 
sions: 



I 

Another related example of this phenomenon (without a 
hard gap) is the problem of a strong Coulomb impurity 
in monolayer graphene that has acquired much interest 
recently..^^'22 
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FIG. 40: [color online] Plot of Im[_D2 (t^)] as a function of the 
frequency ui inside of the band gap for m = 0, r = 10 a, 
V — 50meV and 71 = 72 = 7. Left: from top to bottom 
7 — 1.9, 2.0, 2.1, 2.2; Right: zoom in near the lower band 
edge, from top to bottom 7 = 2.1, 2.3, 2.5, 2.7. 



The important case of a screened Coulomb potential 
generally requires a different approach. Nevertheless, we 
do not anticipate any qualitative discrepancies between 
a potential well and a screened Coulomb potential. We 
expect the screening wave vector to be roughly propor- 
tional to the density of states at the Fermi energy; and 
once the range and the strength of the potential have 
been estimated a potential well can be used to approxi- 
mate the binding energies. We also note that the asymp- 
totic behavior in Eq. (j78p is quite general for a decaying 
potential. 



= Bet[^K,m{uJ,K+,R), ^K,m{(^,K^,R), 'I'/^m (w, /5+ , i?) , 5'/^,„(a;,K_,i?)], 

= Det[^'K,m(w,K+,-R), ^'K,m(w,K-,^), * J,m (t^ , P+ , ^) , (<i, p_ , i?)] , 

= Det[^'K,m(t^,K+,-R), ^'if,m(cJ,K_,i?), * J,.^ (ti, p+ , i?) , ^'./^m (cD, p_ , i?)] , 
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A. Polarization function 

Since we were just discussing the issue of screening 
it fits well to briefly discuss the issue of the dielectric 
function in the biased bilayer (see also Ref. [zl). With 
the introducing of the symmetric (+) and antisymmetric 
(— ) densities (see e.g. Ref. 



(1 



p±(q) = ^^^(k + q) 



\ 
10 
±1 

Vo ±1/ 



7/^(k). (108) 



The usual manipulations then gives the retarded response 
in the symmetric channel as^^ 

X++i~q,u;)^2j2j ^k/(k)«p(k + q)P 

Eiik)-Ei,(k + q)+u + iS- ^ ' 

Here vi (k) is the spinor wave function of band I at mo- 
mentum k. At half filling this expression only have con- 
tributions from I I' and since the wave function overlap 
at q = between different bands for a fixed value of k 
is zero we conclude that x++(qjO) ^ i^i the limit of 
q — > 0. The expression is expected to be dominated by 
the transitions between the E'-.- and £'+^_-bands lead- 
ing to X++(q5 0) — 2g^/F. The Random Phase Ap- 
proximation (RPA) dielectric function is then given by 
e(q) « 1 — (27re^/g)x++(q, 0) which imply that e(q) ~ 1 
as g ^ 0. From this we can conclude that the BGB is 
unable to contribute to the screening of the long-range 
part of the Coulomb interaction. Note that the dimen- 
sionality of the system is crucial for this argument. In 
three dimensions, where the Coulomb interaction goes as 
l/q^ the same argument as in the above usually gives 
a large contribution to e for a semiconductor.^'^ For the 
unbiased bilayer at fi = in the low-energy approxima- 
tion of Eq. ([9]) one finds (using RPA) a screening wave 
vector that is proportional to t^."^^ This is in agreement 
with what one expects for an electron gas in 2D where 
the screening wave vector is proportional to the the effec- 
tive mass. For a more detailed discussion of the unbiased 
graphene bilayer dielectric function including the trigonal 
warping see Ref. [tI. 



COHERENT POTENTIAL 
APPROXIMATION 



XV. 



As discussed above in Sec. lXIIl for a finite density ni of 
impurities, the bound states can interact with each other 
leading to the possibility of band gap renormalization and 
the formation of impurity bands. A simple, but crude, 
theory of these effects is the CPA.~i22, In this approxima- 
tion, the disorder is treated as a self-consistent medium 



with recovered translational invariance. The medium is 
described by a set of four local self-energies which are 
allowed to take on different values on all of the inequiv- 
alent lattice sites in the problem. In fact this section 
is a straightforward extension to the biased case of the 
methods applied in Sec. |V] for the unbiased case. The 
self-energies are chosen so that there is no scattering on 
average in the effective medium. It has been argued that 
the CPA is the best single-site approximation to the full 
solution of the problem.— 

In the following we often suppress the frequency de- 
pendence of the self-energies for brevity. The expression 
for the diagonal elements of the Green's function G is 
given in App. [d1 We follow the standard approach to 
derive the CPA (see for example Refs. [sHItsI I). and we 
obtain the self-consistent equations: 



1 ^ {U — Sqj)Gqj 



(110) 



The limit of site dilution (or vacancies) used in Sec. |V] 
is obtained in the limit U oo leading to the self- 
consistent equations: 



rij 

Gaj 



(111) 



An explicit expression for the local propagators Gaj is 
given in App. |DJ Using the expressions obtained there 
the self-consistent equations for U oo becomes: 



-Gai 



— = -Gbi 

^Bl 



/3i (a -"2/32^0), (112a) 



From these equations it is straightforward to obtain the 
density of states (DOS) on the different sublattices aj 
from Paji^j) — — ImGctj((jj + id)/TT. In the clean case, 
one finds: 



LO - V/2 



2A2 



LoVi2-x) 



:_3a) 



, ti{u: + V/2){2-x) 
Pai 



(113b) 



for \uj\ > Eg/2. Here x = (0,1,2) for (|w| < V/2, 
V/2 < |w| < ^tl + Vy4, y/tl + Vyl < |w|). The 
corresponding quantities in plane 2 are obtained by the 
substitution V ^ —V. In the limit of y ^ we re- 
cover the known unbiased result of Eq. . Notice that 
the square-root singularity starts to appear already above 
V/2 on the Bl sublattice. There is also a divergence on 
the Al sublattice but the coefficient in front is usually 
much smaller. The DOS on the Al sublattice vanishes 
at a; = V/2 while the DOS on the Bl sublattice is finite 
there. 

The numerically calculated density of states for U 
oo is shown in Fig. 1411 The impurity band evolves from 
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FIG. 41: [color online] Left: DOS as a function of the energy 
(in units of Eg) close to the conduction band edge for differ- 
ent impurity concentrations (see inset), U = —1 eV. Right: 
Details of the DOS inside of the gap for different impurity 
concentrations for ?7 ^ oo. In both cases V = 40meV. 



the single-impurity B2 bound state which for the param- 
eters involved is located at e w 0.3i?g. Further evidence 
for this interpretation is that the total integrated DOS 
inside the split-ofF bands for the two lowest impurity con- 
centrations is equal to Ui. It is worth mentioning that 
the width of the impurity band in the CPA is likely to be 
overestimated. The reason for this is that that the use 
of effective atoms, all of which have some impurity char- 
acter, increases the interaction between the impurities.2S 
For smaller values of the impurity strength the single- 
impurity bound states are all weakly bound (cf. Fig. [37|) 
and the "impurity bands" merge with the bulk bands as 
shown in Figs. |4T] and |42l The bands have been shifted 
rigidly by the amount riiV for a more transparent com- 
parison between the different cases. The smoothening 
of the singularity and the band gap renormalization is 
apparent. Observe also that the band edge moves fur- 
ther into the gap at the side where the bound states 
are located. It is likely that the CPA gives a better ap- 
proximation for these states since by Eq. (|75|) they are 
weakly damped almost propagating modes. Notice that 
the gap and the whole structure of the DOS in the re- 
gion of the gap is changing with V ^ and in particular 
the possibility that the actual gap closes before y = 
because of impurity-induced states inside of the gap. Fi- 
nally we note that this observation is consistent with the 
results of numerically exact calculations using the recur- 
sive Green's function method for strong disorder iS 



XVI. EFFECTS OF TRIGONAL DISTORTIONS 

Before we conclude this paper we would like to briefly 
comment on the effect of the 73-term on our findings in 
the previous sections. The effects of 73 on the spectrum 
of the BGB was discussed in Sec. lIIIl where it was shown 
that this term breaks the cylindrical symmetry and leads 
to the "trigonal distortion" of the bands. In the BGB 
the result is three copies of a more conventional ellip- 
tic dispersion at the lowest energies near the band edge. 
Using the same method as in Section IXII we find that. 



FIG. 42: [color online] DOS as a function of energy (in units 
of Eg) for different values of the applied bias V (see inset) 
and [/ = -10 eV. Left: m = 10"^; Right: n,; = 10"^ 



also for an elliptical band edge, a Dirac delta potential 
always generates a bound state in 2D. The divergence of 

G is generally only logarithmic as the band edge is ap- 
proached however, whereas the divergence is an inverse 
square root without 73. More confined bound states with 
larger binding energies sample a larger area of the BZ. 
Therefore we do not expect that the small details at the 
band edge to significantly modify the results that we ob- 
tain with the minimal band model for these states. 

Another observation is that when there is a finite den- 
sity of impurities in the BGB the self-energies can become 
quite large as we have seen in the previous section. Con- 
sider the case T/ = 50 meV, for which iy-E^ /E^ « 0.01. 
Therefore, by looking at Fig. O we see that 73 smooth 
out the square-root singularity on this scale. Compar- 
ing with Fig. m] we see that for an impurity of strength 
U = —1 eV the trigonal distortion would correspond to 
a density of impurities of around nj ^ .001. In the case 
that the gap is filled up with impurity-induced states (see 
Fig. l42)) . the disorder-induced energy scale is much larger 
than that generated by 73. Therefore we argue that the 
possibility that the gap closes before y = is robust to 
the presence of a 73, even if it is as large as the values 
quoted in the graphite literature. 



XVII. CONCLUSION AND OUTLOOK 

Graphene research is one of the fastest growing fields 
in condensed matter research since the isolation by 
Novoselov et al. of the first graphene flake in 20041^ 
Three years after that, and after hundreds of theoret- 
ical papers on the subject,"^ the physics of single layer 
graphene is relatively well understood and very few con- 
troversies remain. Meanwhile, the study of multilayer 
graphene, and particularly bilayer graphene, continues 
to be, experimentally and theoretically, an open field of 
research. Partially, this can be assigned to the natu- 
ral attraction of researchers to the "one atom thin" ma- 
terial and its unusual electronic spectra. Nevertheless, 
graphene bilayer is equally thin (two atoms thick, indeed) 
and has also 2D Dirac spectrum (albeit massive) with un- 
usual properties. Graphene bilayer is also more prone to 
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show strong electron-electron correlation effects, such as 
magnetism^ and charge density waveSfl phases, because 
of its finite density of states at the Dirac point, unlike its 
single layer counterpart where interactions are at most 
marginal in renormalization group sensei^ Interestingly 
enough, the study of the effect of electron-electron inter- 
actions in graphene is a field in its infancy i'^^'^^ 

Furthermore, graphene bilayer is equally easy (or hard) 
to find or produce epitaxially. Besides its intriguing 
electronic properties, the graphene bilayer is a promis- 
ing candidate to bulk electronic devices with proper- 
ties that are insensitive to surface (edge) defects such 
as graphene nanoribbons and quantum dots. Perhaps 
even more interesting is the fact that graphene bilayer 
is the only known material that has an electronic gap 
between conduction and valence bands that can be fully 
controlled by the application of a transverse electric field 
(a tunable gap semiconductor), as has been demonstrated 
experimentally,^^"^'^ This property opens up an enormous 
number of possible ways to use bilayer graphene, from 
transistors to lasers working in the terahertz regime. 

Nevertheless, in order to be able to use graphene bi- 
layer (and multilayers) in device applications, one has 
to understand how material issues, such as disorder, af- 
fect its electronic properties. This was the main aim of 
this work, namely, to understand how disorder affects the 
electronic properties of bilayer (and multilayer) graphene 
in its most basic model. We have shown that the elec- 
tronic self-energies can be calculated analytically within 
the T-matrix and CPA presenting some unusual features 
that can be measured either by transport or spectroscopy 
(ARPES and STM). We have calculated a series of im- 
portant physical properties such as spectral functions, 
and frequency dependent conductivities. We have also 
studied the problem of bound states in the biased bilayer 
graphene and their effect in the electronic structure and 
have shown that the properties of these bound states can 
be equally well controlled by applied transverse fields. 
We also point that we have left open issues associated 
with trigonal distortions. At this point in time, it is not 
clear that such effects, associated with 73, are going to 
be the same as observed in 3D graphite and more ex- 
perimental studies are needed in order to investigate the 
issue. We hope that our results will stimulate more ex- 
perimental studies of these amazing new materials. 
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APPENDIX A: MINIMAL CONDUCTIVITY OF 
THE GRAPHENE BILAYER INCLUDING 
TRIGONAL WARPING 

In this appendix we provide an alternative derivation 
of the value of the minimal conductivity of the graphene 
bilayer in the presence of trigonal warping. The conduc- 
tance of a wide strip of graphene at the Dirac point is 
mediated by evanescent modes that connect the leads. 
We define the hamiltonian as 



n = 





vpiK - iky) 



vpikx + iky) 




(Al) 



and use the Landauer formalism described in Refs. l83ll84l . 
We take the width W of the sample to be much larger 
than its length L. If we assume that the leads on the 
right and on the left are heavily doped clean graphene, 
the incoming, reflected, and outgoing waves can be ap- 
proximated as: 



ref 




ikyy 



ikyy 



ikyV 



(A2a) 
(A2b) 
(A2c) 



where t{ky) [r{ky)] is the transmission [reflection] ampli- 
tude of the mode with perpendicular momentum ky. The 
wavefunction in the central region, < a; < L, at zero 
energy, can be written as 



^1 = A 



iky y 







B 







JkyV 



(A3) 



The matching conditions at the edges at x = and x — L 
are: 



r{ky) 
r{ky) 

t{ky) 
t{ky) 



= \/2A, 
= \/2B, 
= V2Ae 



resulting in the transmission probability 

T{ky) = \t{ky)\^ = ^ . 

cosh {kyL) 

The conductance per channel is thus given by 



^ ~ 1i2tt ' '^^vT{ky) 



(A4a) 
(A4b) 
(A4c) 
(A4d) 

(A5) 
(A6) 



We will now extend this result to the anisotropic Dirac 
equation. The Hamiltonian is 







fJx kx 



Vxkk 



1' ky 



(A7) 
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where the Fermi velocities {vx and Vy) are allowed to be 
different in the x and y directions. We use the incoming 
and outgoing wavefunctions wavefunctions in Eq. (jA2|l . 
and generalize the wavefunctions in the graphene junc- 
tion, Eq. ([X3|) to 



^ = A 



JkyV 







B 







JkyV 



(A8) 



so that 



G 



-)(-) 



W 

T' 



(AlO) 



where the Dirac equation implies that 

Matching the wave functions at the contacts with the 

leads, we find that the generalization of Eq. (jASp is 



T{ky) = 



1 



cosh" 



( ^y ky 



(A9) 



If the junction is rotated by an angle 9 with respect to the 
main axes of the anisotropic Dirac equation, the Hamil- 
tonian becomes: 



n 





Va[kx cos{9) — ky sm{6)] — ivh[kx sin(6') + ky cos(6')] 



Va[kx cos{9) — ky sin(6')] + ivij[kx sin(6') + ky cos( 




iAll) 



where Va and Vb are the Fermi velocities along the two 
principal axes. The wave function in the central region 
is now 



^ = A 



where 



VaVb 







k' 



vlcos^9)+vlsm^{9) 
sm{9)cos{9){vl-vl) ^ 



vl cos'^{9) 
Using this we also obtain 

T{ky 



■sm^{9) 



1 



cosh 



VaVfykyL 

Va cos^{e)+v'^ sin^(e) 



leading to 



G 



vl cos^{9)+vl siT?{9)\W 



VaVb 



JkyV 

(A12) 
(A13a) 
(AlSb) 

(A14) 
(A15) 



In a graphene bilayer, including trigonal warping but 
ignoring terms that couple sites in the same sublattice, 
we have four Dirac points. One of them is isotropic, 
with Va = Vb, and the three others are anisotropic, 
with Vb — 3va- The principal axes at these three Dirac 
points form angles with respect to a barrier which can be 
parametrized as 0o , + 27r/3 and 0o + 47r/3 , where de- 
pends on the orientation of the barrier. The conductance 
is therefore given by 



2 \Vb Va 



(A16) 



This expression is independent of the angle ^o- For 
Vb/va = 3, we find G = 6 x [e^/ (7rfe)] x (W/L) per channel, 



in agreement with Refs. 



APPENDIX B: DENSITY OF STATES IN 
MULTILAYER GRAPHENE 

In this appendix we derive explicit expressions for the 
DOS in graphene multilayers. The expressions are used 
in Section IVll To calculate the DOS in graphite we must 
perform two integrals to get G. One integral we need is 



h = 



1 



1 



2D. 



2D4 



(Bl) 



First we perform the perpendicular integral using com- 
plex variables to rewrite the integral as a contour integral 
around the unit circle and then picking up the pole inside: 



7r/2 



dk I 



t/2 



1 1 

1 



dz 



i tj_{z'^ + 1)- Az 



_ 2tt sign[Re(A)] 

where A = p^/ (tj — Ss) — (tj — Sa). Note that the func- 
tion sign[Re(^)] changes sign just where the branch of 
the square root does. Moreover the square root is purely 
imaginary there. Therefore the function is actually con- 
tinuous across the point where Re(A) = 0. From now on 
in this appendix, as well as in Appendix IC 2i we choose 
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the branch of the square root such that the sign of the real 
part is included, with this convention A and — 4^^ 
always lies in the same quadrant of the complex plane. 
Because of the form of A we can use the integral formula 



e\v')= dip') 



1 



= (cj - SB)log 



A 



- 



, (B2) 



directly, since the argument of the log does not cross any 
branch cut. Thus the result of the integral is 



log 



A(A2) + y/A^{A^)-4tl 



A{0) + ^A^0)-4t 
Finally to leading order in A we get 

2A2 



(B3) 



Ji - log 
Similarly the integral 

h - 



2 

(B4) 



2n dfc [-2i_LCOs(fc) , 2tj_cos{k) 

a[p ) 



10 J-n/2 

can be written as 
1 



2D. 



A^ 



A" 



A 



2 -4l^ 



Now we may use the integral formula 
A 



dip' 



iuj 



A'-4tl), 



(B5) 
(B6) 

(B7) 



and one can once again convince oneself that there are no 
contribution from the possible crossing of the branch at 



J 



Re(A) = 0. With the help of this the resulting expression 
becomes 



1 



UJ — 'Si 



-{a2 -iiu- Ss) [^/A^il^)^t 



A2(0)-44]}. (B8) 



Finally, keeping only the leading term in the expansion 
in A we get 



/2 = (cj - Y.a) 



ic 



(B9) 



APPENDIX C: CONDUCTIVITY KERNELS 

In this appendix we derive formulas for the conductiv- 
ity kernels that we use in Sec. I Villi and Sec. lIXI First we 
rewrite the kernels with the identity 



Im[Gi(e)] Im[G2(e + cj)] 

= i Re [Gtie)Gfie + u) - Gf (e)G^ (e + u:)] 



iRe[^ 7Gi(Si - i^T i)G2iE2 + 1T2 



7=± 

and introduce the notations 



, (CI) 



Ea ^ e-- SA(e), 

Eb = e - SB(e), 

Ea = e + tt) - SA(e + t^), 

Eb = e + — SB(e + t^)- 



(C2a) 
(C2b) 
(C2c) 
(C2d) 



1. Bilayer graphene 



The integrals we need for the bilayer are easy to obtain since there is never any problems with branches of the 
logarithms. The kernels can be expressed in terms of the following integrals: 



1 



EBiEA+f3h 



'^t:^ EBiEA + at^)-EBiEA+f3t±) l-EeiEA + at^) 



■log 



~EBiEA + l3t^) 



(C3) 



1 



EeiEA + at^) 



4 ^ EBiEA + fiti_) - Eb{Ea + atA_) 



log 



^EBiEA + ati_) 
-EBiEA + f3t^) 



(C4) 



32 



A2 



-SB(^A+/3ij.)l0g 



A2 



}, (C5) 



and 



db^) [g2A(e)g?A(e + ^) - gA2(^)gAS(^ + ^)] 



EbEb 



log 



2 ^ Eb{Ea + - £;B(i;A - + atj^) 



-EB{EA~at^) 



(C6) 



In fact, although the cutoff A enter the expression in Eq. (|C5p . the final result is actually independent of A. For the 
DC conductivity it is convenient to work out that for two retarded propagators we have the relation 



rf(p^)[gSAWggB(e)+g^g(e)g2gW] 



(C7) 



2. Multilayer graphene 

For multilayer graphene we have to perform two integrals, they are 

r-/2 dk ~ - 







7r/2 



D „D _L „D D I 9o,ND ND] 
AAgBB + gBBgAA + ^gABgABj 



(C8) 



and 



dk 



-7r/2 



gMA+gMS] Bin^(fc)- 



(C9) 



Exactly like in the case above we find Ji = —2 when w = and both the propagators are retarded. First we perform 
the perpendicular integral using a contour integral as we did when we computed the DOS in Appendix [B) In the 
following the branch of the square root includes the sign of the real part ^ — 4t^ = sign[Re(y4)]-y/A2 _ 4^2^^ This 

implies that A and ^ — \t\ always lies in the same quadrant. The results we need are 



dk_ 

-7r/2 7^ 



gAA(e)gEB(e + ^) 



2Eb I _ 



V 

Eb 



A~A A+A 



^A^ - 4t 



P 

Eb 



A-A A+A 



- itl 



(CIO) 



'"'/^ dk 

-g^g(e)g^g(e + - 

-7T/2 



2EbEb ^ 



l- 1 ^ 1 ] 


1 


I 1 ^ 1 1 


U-I A + Ai 


^A^-iti^ 


^A-A A + Ai 



— |, (cii) 

< ,9 J 



- 44 



and 



JT/2 V [gAA(e)g?A(^ + -) + gXS(e)g2S(e + c.)] sin^(fc) = -L{_i + _1 



>1 



^2 _ 4^2^ 



yl2 - 44 



Where 



^ = pV^s - Ea, 
A = p^Eb-Ea. 



(C12) 



(C13) 
(C14) 
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Adding up all the contributions for the Ji we get after some rearrangements 



+ 



2rf_ 



2d+ 



1 + 



1 + 



K - c+ 



1 1 

+ 



— c_ 



(C15) 



where 



d± 



= hiBJiiB— 



d± 



En dz Ef 



I ^ I _ Eb ±Eb 
Eb Eb Eq Eq 



(C16) 
(C17) 



It is convenient to define the integral 

1 



I ^ |log(p' - c±) - log[Ai?± + JA^-AtlJBl-Atl-Ul] |, (C18) 



in which we have introduced 



B± = c±/Eb — Ea, 



(C19) 



and the square roots are again chosen such that B± and ^ B^. — are in the same quadrant in the complex plane. 
Using this we may write 

/ Ja^ - Atl = Ja^ - Atl + B± log [A + jA^-4tl] + [Bi - 4ti] {p% (C20) 

If we use this formulas blindly and just plug in the endpoints the resulting expressions are 



Ji = 



d^_d^ 1_ 

2d^~ 2dZ' Wn 



+ 



1 



2d+ 2d- Eb 
d?_c+ 



+ 



2d+ 



die. 



2d- 



(C21) 



and 



= M-^^T. [/^- + i^^^^H.) - ^/^\z) 



[Bi-Ati]e\z)-[Bi-Ati]e\z) 



2U(2), 



(C22) 



z=0 



One must be careful with the imaginary part of ^^^^ how- titles to write: 
ever. First wc note that _B+ = — i?+ and B^ = _B_, 
which implies that the log(p^ — c±) term in ^(^^ does jgg 
not contribute. Secondly, we write A = 2t± cosh(a) and 
B = 2t± cosh(/3) and use hyperbolic trigonometric iden- 



AB + ^ A^ - At\^ B"^ - At]_ - Atl 
= log{4ti[cosh(a + /3) - 1]} 



log(8ti) + 21og[sinh(^i^)] . (C23) 
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By convention Re(a) > and — tt < Im(a) < tt (and 
the same goes for /3). Therefore — tt < Ini(Q; + /3)/2 < tt 
and the argument of the log never crosses the branch 
cut along the negative real axis. Therefore the repre- 
sentation of A and B in terms of hyperbolic functions 
automatically takes care of the phase information of the 
argument. Finally, the case when B is purely real and 
negative (this is relevant for the calculation of the DC 
conductivity) requires that one should choose Im(/3) < 
so that sign[Re(B)] y^B'^ - U\ = 2tj_ sinh(/3). 



APPENDIX D: LOCAL PROPAGATOR IN THE 
BIASED GRAPHENE BILAYER 

In this appendix we provide some details on the calcu- 
lation of the local propagator in the BGB that we use in 
Sec. IXVI Introducing the notations: 



where we have defined 



z± = ± o V ("iPi - a2P2)^ + 4tj_/3i/32. 



2 2 
Introducing the integrals 



- In 



In 



A2 



and 



A2- 



In 



A2 



(D4) 



A_2 M 



ai 


= LU 


- V/2- 




(Dla) 


Pi 


= LU 


- V/2- 




(Dlb) 


a2 


= LU 


+ V/2- 




(Die) 


P2 


= LU 


+ y/2- 


Sb2, 


(Did) 



one can write 

Gaiai 
Gbibi 



/3i(a2/32 



D 

ai{a2/32 - k^) - tl/32 
D 



(D2a) 
(D2b) 



The equations for the corresponding quantities in layer 
2 are obtained by exchanging the indeces 1 2 every- 
where. The denominator can be written as 



D = (ai/3i-fc")(a2/32-fc^)-il/3i/32 = (r-z_)(r-z+), 

(D3) 



In 



A2-Z+ 



In 



A2 



— z_ In 



z_ In 



A2_ 

— z_ 



, (D6) 



we can easily compute G. Using the explicit form of 
the propagators in Eq. (|D2p and the same continuum 
approximation as in Eq. (|25p we obtain 

-Gai = /3i (6 -"2/32^0), (D7a) 
-Gbi - ai{^i~a2f32^Q)+tlMo- (D7b) 

Again the corresponding quantities in plane 2 are ob- 
tained by exchanging the indeces 1 <-!■ 2 everywhere. 
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